
#include <avtMDSCHISMFileFormatImpl.h>

#include <string>
#include <iostream>
#include <sstream>
#include <time.h>
#include <math.h>

#include <vtkCharArray.h>
#include <vtkDoubleArray.h>
#include <vtkFloatArray.h>
#include <vtkIntArray.h>
#include <vtkLongArray.h>
#include <vtkShortArray.h>
#include <vtkUnsignedCharArray.h>
#include <vtkUnsignedIntArray.h>
#include <vtkUnsignedShortArray.h>

#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkCellType.h>
#include <vtkFieldData.h>
#include <vtkInformation.h>
#include <vtkRectilinearGrid.h>
#include <vtkStreamingDemandDrivenPipeline.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>

#include <avtDatabaseMetaData.h>
#include <avtVariableCache.h>
#include <avtGhostData.h>
#include <DBOptionsAttributes.h>
#include <Expression.h>

#include <InvalidVariableException.h>
#include <InvalidDBTypeException.h>
#include <InvalidTimeStepException.h>
#include <InvalidFilesException.h>
#include <DBYieldedNoDataException.h>
#include <DebugStream.h>
//L3 #include <malloc.h>
#if defined(__MACH__)
#include <stdlib.h>
#else
#include <malloc.h>
#endif

#include <BufferConnection.h>
#ifdef PARALLEL
#include <avtParallel.h>
#endif

#include "MDSCHISMMeshProvider.h"
#include "SCHISMFileUtil10.h"
#include "Average.h"
#include "MeshConstants10.h"
#include "MDSCHISMOutput.h"
#include "avtMDSCHISMFileFormat.h"
#include "Registar.h"

using     std::string;
using     std::stringstream;

const std::string NODE = MeshConstants10::NODE;
const std::string FACE = MeshConstants10::ELEM;
const std::string SIDE = MeshConstants10::EDGE;
const std::string UNKOWN = "unkown";
const int NODESPERELE = MeshConstants10::MAX_NUM_NODE_PER_CELL;
const int NODESPERWEDGE = NODESPERELE * 2;


avtMDSCHISMFileFormatImpl::avtMDSCHISMFileFormatImpl() :
	m_initialized(false),
	m_mesh_is_static(true),
	m_plugin_name("SCHISM_output_plugin"),
	m_data_description("data_description"),
	m_mesh_var("Mesh2"),
	m_var_label_att("long_name"),
	m_var_location_att("location"),
	m_mesh_3d("3D_Mesh"),
	m_layer_mesh("Layer_Mesh"),
	m_mesh_2d("2D_Mesh"),
	m_mesh_2d_no_wet_dry("2D_Mesh_no_wet_dry"),
	m_side_center_point_3d_mesh("side_center_3D"),
	m_side_center_point_2d_mesh("side_center_2D"),
	m_face_center_point_3d_mesh("face_center_3D"),
	m_dim_time(MeshConstants10::DIM_TIME),
	m_time(MeshConstants10::TIME),
	m_node_depth(MeshConstants10::NODE_DEPTH),
	m_node_depth_label(MeshConstants10::NODE_DEPTH),
	m_dim_layers(MeshConstants10::DIM_LAYERS),
	m_dim_var_component(MeshConstants10::DIM_VAR_COMPONENT),
	m_time_ptr(NULL),
	m_node_x_ptr(NULL),
	m_node_y_ptr(NULL),
	m_node_z_ptr(NULL),
	m_kbp_data(NULL),
	m_kbp_node_filled(false),
	m_kbp_side_filled(false),
	m_kbp_ele_filled(false),
	m_cache_kbp_id(-1),
	m_surface_state_suffix("_surface"),
	m_bottom_state_suffix("_near_bottom"),
	m_depth_average_suffix("_depth_average"),
	m_dry_surface(MeshConstants10::DRY_SURFACE),
	m_dry_wet_flag(0)
{
	// AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
	m_center_map[NODE] = AVT_NODECENT;
	m_center_map[FACE] = AVT_ZONECENT;
	m_center_map[UNKOWN] = AVT_UNKNOWN_CENT;
	m_var_name_label_map[m_node_depth_label] = m_node_depth;
	
}

FileFormatFavorInterface * avtMDSCHISMFileFormatImpl::create()
{
	return new avtMDSCHISMFileFormatImpl();
}


// ****************************************************************************
//  Method: avtEMSTDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Wed Mar 13 09:13:49 PDT 2013
//
// ****************************************************************************

int
avtMDSCHISMFileFormatImpl::GetNTimesteps(const std::string& a_filename)
{
	Initialize(a_filename);
	//debug1 << "num time step is"<<m_num_time_step <<"\n";
	return m_num_time_step;
}

//void  avtMDSCHISMFileFormatImpl::ActivateTimestep(const std::string& a_filename)
//{
	//Initialize(a_filename);
//}

avtMDSCHISMFileFormatImpl::~avtMDSCHISMFileFormatImpl() 
{ 
	if (m_time_ptr)
	{
		delete m_time_ptr;
	}

	if (m_node_x_ptr)
	{
		delete m_node_x_ptr;
	}
	if (m_node_y_ptr)
	{
		delete m_node_y_ptr;
	}


	if (m_kbp_data)
	{
		delete m_kbp_data;
	}

	std::map<int, int*>::iterator it;
		for (it = m_kbp_node.begin(); it != m_kbp_node.end(); it++)
		{
			if (it->second)
			{
				delete it->second;
			}
		}
		for (it = m_kbp_ele.begin(); it != m_kbp_ele.end(); it++)
		{
			if (it->second)
			{
				delete it->second;
			}
		}
		for (it = m_kbp_side.begin(); it != m_kbp_side.end(); it++)
		{
			if (it->second)
			{
				delete it->second;
			}
		}
		for (it = m_kbp_prism.begin(); it != m_kbp_prism.end(); it++)
		{
			if (it->second)
			{
				delete it->second;
			}
		}

		for (it = m_node_dry_wet.begin(); it != m_node_dry_wet.end(); it++)
		{
			if (it->second)
			{
				delete it->second;
			}
		}

		for (it = m_ele_dry_wet.begin(); it != m_ele_dry_wet.end(); it++)
		{
			if (it->second)
			{
				delete it->second;
			}
		}

		for (it = m_side_dry_wet.begin(); it != m_side_dry_wet.end(); it++)
		{
			if (it->second)
			{
				delete it->second;
			}
		}

		std::map<int, MDSchismOutput*>::iterator it1;
		for (it1 = m_data_files.begin(); it1 != m_data_files.end(); it1++)
		{
			if (it1->second)
			{
				delete it1->second;
			}
		}

		std::map<int, MDSCHISMMeshProvider*>::iterator it2;
		for (it2 = m_external_mesh_providers.begin(); it2 != m_external_mesh_providers.end(); it2++)
		{
			if (it2->second)
			{
				delete it2->second;
			}
		}
}
// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Wed Mar 13 09:13:49 PDT 2013
//
// ****************************************************************************

void
avtMDSCHISMFileFormatImpl::FreeUpResources(void)
{
	return;
}


void avtMDSCHISMFileFormatImpl::set_var_name_label_map(const std::map<std::string, std::string> a_map)
{
	m_var_name_label_map = a_map;
}
void avtMDSCHISMFileFormatImpl::set_var_mesh_map(const std::map<std::string, std::string> a_map)
{
	m_var_mesh_map = a_map;
}
void avtMDSCHISMFileFormatImpl::set_var_horizontal_center_map(const std::map<std::string, std::string> a_map)
{
	m_var_horizontal_center_map = a_map;
}
void avtMDSCHISMFileFormatImpl::set_var_vertical_center_map(const std::map<std::string, std::string> a_map)
{
	m_var_vertical_center_map = a_map;
}
//broadcast string map from rank 0 to all other ranks
void  avtMDSCHISMFileFormatImpl::broadCastStringMap(std::map<std::string, std::string>& a_map, int myrank)
{
#ifdef PARALLEL
	std::vector<std::string> key_list;
	std::vector<std::string> val_list;
	if (myrank == 0)
	{
		for (std::map<std::string, std::string>::iterator it = a_map.begin(); it != a_map.end(); ++it)
		{
			key_list.push_back(it->first);
			val_list.push_back(it->second);
		}
		BroadcastStringVector(key_list,myrank);
		BroadcastStringVector(val_list,myrank);
	}

	if (myrank != 0)
	{
		std::vector<std::string>::iterator key_it = key_list.begin();
		std::vector<std::string>::iterator val_it = val_list.begin();
		while ((key_it != key_list.end()) && (val_it != val_list.end()))
		{
			a_map[*key_it] = *val_it;
			key_it++;
			val_it++;
		}
	}
#endif
}


void avtMDSCHISMFileFormatImpl::BroadcastGlobalInfo(avtDatabaseMetaData *metadata)
{
#ifdef PARALLEL
	int rank = PAR_Rank();
	
	//MPI_Comm_rank(VISIT_MPI_COMM, &rank);
	//
	// Broadcast Full DatabaseMetaData
	//
	BufferConnection tmp;
	int n;
	if (rank == 0)
	{
		metadata->SelectAll();
		n = metadata->CalculateMessageSize(tmp);
		debug1 <<" n is "<<n<<"\n";
		
	}
	BroadcastInt(n);
	debug1 << "after update n is " << n << "\n";
	unsigned char *buff = new unsigned char[n];
    int *buffint = new int[n];
	debug1 << "mpi bcast \n";
	if (rank == 0)
	{
		
		metadata->SelectAll();
		metadata->Write(tmp);
		
		for (int i = 0; i < n; i++)
		{
			tmp.Read(&buff[i]);
			buffint[i] = buff[i];
			//debug1 <<" "<<buff[i];
		}
		 //debug1 << "\n";
		
	}
	debug1 << "mpi bcast1 \n";
	BroadcastIntArray(buffint, n);
	//MPI_Init(NULL, NULL);
	//MPI_Bcast(buff, n, MPI_UNSIGNED_CHAR, 0, VISIT_MPI_COMM);
	debug1 << "mpi bcast2 \n";
	if (rank != 0)
	{
		for (int i = 0; i < n; i++)
	   {
			buff[i] = buffint[i];
			//debug1 << " " << buff[i];
		}
		//debug1 << "\n";
		tmp.Append(buff, n);
		metadata->Read(tmp);
	}
	delete[] buff;
	delete[] buffint;
	broadCastStringMap(m_var_name_label_map, rank);
	broadCastStringMap(m_var_mesh_map, rank);
	//MPI_Finalize();
#endif
}

// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Wed Mar 13 09:13:49 PDT 2013
//
// ****************************************************************************

void
avtMDSCHISMFileFormatImpl::PopulateDatabaseMetaData(avtDatabaseMetaData *a_metaData, const std::string a_data_file, int a_timeState)
{
	//debug1 << "begin fill metadata";
#ifdef PARALLEL
	if (PAR_Rank() == 0)
	{
		debug1 << "not root no meta filling job on rank"<<PAR_Rank()<<"\n";
		
		Initialize(a_data_file);
		//
		avtCentering  nodeCent = AVT_NODECENT;
		avtCentering  zoneCent = AVT_ZONECENT;
		string mesh_name = m_mesh_3d;

		// AVT_RECTILINEAR_MESH, AVT_CURVILINEAR_MESH, AVT_UNSTRUCTURED_MESH,
		// AVT_POINT_MESH, AVT_SURFACE_MESH, AVT_UNKNOWN_MESH
		avtMeshType mt = AVT_UNSTRUCTURED_MESH;
		//
		int nblocks = m_number_domain;
		//debug1 << "num of domain "<<m_number_domain<<"\n";
		int block_origin = 0;
		int cell_origin = 0;
		int group_origin = 0;
		int spatial_dimension = 3;
		int topological_dimension = 3;
		double *extents = NULL;
		int * bounds = NULL;
		//add node center 3d mesh
		if (m_external_mesh_providers.begin()->second->provide3DMesh())
		{
			avtMeshMetaData *mmd = new avtMeshMetaData(bounds,extents, mesh_name,
				nblocks, block_origin, cell_origin,group_origin,spatial_dimension,topological_dimension, mt);
			a_metaData->Add(mmd);
		}

		//
	   // add layered 2d mesh
		mesh_name = m_layer_mesh;
		topological_dimension = 2;

		if (m_external_mesh_providers.begin()->second->provide3DMesh())
		{
			avtMeshMetaData *mmd = new avtMeshMetaData(bounds, extents, mesh_name,
				nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
			a_metaData->Add(mmd);
		}

		//
		// add surface 2d mesh
		//
		mesh_name = m_mesh_2d;

		spatial_dimension = 2;
		topological_dimension = 2;

		avtMeshMetaData *mmd1 = new avtMeshMetaData(bounds, extents, mesh_name,
			nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
		a_metaData->Add(mmd1);

		//add 3d side center point mesh
		mesh_name = m_side_center_point_3d_mesh;

		spatial_dimension = 3;
		topological_dimension = 0;

		if (m_external_mesh_providers.begin()->second->provide3DMesh())
		{
			avtMeshMetaData *mmd = new avtMeshMetaData(bounds, extents, mesh_name,
				nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
			a_metaData->Add(mmd);
		}

		//add 3d  face center point mesh
		mesh_name = m_face_center_point_3d_mesh;

		spatial_dimension = 3;
		topological_dimension = 0;

		if (m_external_mesh_providers.begin()->second->provide3DMesh())
		{
			avtMeshMetaData *mmd = new avtMeshMetaData(bounds, extents, mesh_name,
				nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
			a_metaData->Add(mmd);
		}

		// add 2d side center point mesh
		mesh_name = m_side_center_point_2d_mesh;

		spatial_dimension = 2;
		topological_dimension = 0;

		avtMeshMetaData *mmd2 = new avtMeshMetaData(bounds, extents, mesh_name,
			nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
		a_metaData->Add(mmd2);

		//
		// add water surface and depth scalar
		string mesh = m_mesh_2d;
		avtScalarMetaData *smd1 = new avtScalarMetaData(m_node_depth_label,mesh, nodeCent);
		a_metaData->Add(smd1);



		// m_var_mesh_map[m_node_surface_label] = mesh;
		m_var_mesh_map[m_node_depth_label] = mesh;

		//add 3D node level label
		string mesh3d = m_mesh_3d;
		avtScalarMetaData *smd2 = new avtScalarMetaData(MeshConstants10::NODE_LEVEL, mesh3d, nodeCent);
		a_metaData->Add(smd2);

		//add 3D element level label
		mesh3d = m_layer_mesh;
		avtScalarMetaData *smd3 = new avtScalarMetaData(MeshConstants10::ELE_LEVEL, mesh3d, zoneCent);
		a_metaData->Add(smd3);

		//add 3D side level label
		mesh3d = m_side_center_point_3d_mesh;
		avtScalarMetaData *smd4 = new avtScalarMetaData(MeshConstants10::SIDE_LEVEL, mesh3d, nodeCent);
		a_metaData->Add(smd4);

		//add 3D layer lable (for prism center data)
		mesh3d = m_mesh_3d;
		avtScalarMetaData *smd5 = new avtScalarMetaData(MeshConstants10::LAYER, mesh3d, zoneCent);
		a_metaData->Add(smd5);


		//add node bottom label
		mesh = m_mesh_2d;
		//a_avtFile->addScalarVarToMetaData(a_metaData,  MeshConstants10::NODE_BOTTOM, mesh, nodeCent);
		//a_avtFile->addScalarVarToMetaData(a_metaData,  MeshConstants10::FACE_BOTTOM, mesh, zoneCent);

	   //a_avtFile->addScalarVarToMetaData(a_metaData,  MeshConstants10::EDGE_BOTTOM, mesh, nodeCent);


		PopulateStateMetaData(a_metaData, a_timeState);
		debug1 << "finish populate metadata \n";
		
	}
	//Barrier();
	BroadcastGlobalInfo(a_metaData);
	debug1 << "finish populate and broadcast metadata \n";
#else
	Initialize(a_data_file);
	//
	avtCentering  nodeCent = AVT_NODECENT;
	avtCentering  zoneCent = AVT_ZONECENT;
	string mesh_name = m_mesh_3d;

	// AVT_RECTILINEAR_MESH, AVT_CURVILINEAR_MESH, AVT_UNSTRUCTURED_MESH,
	// AVT_POINT_MESH, AVT_SURFACE_MESH, AVT_UNKNOWN_MESH
	avtMeshType mt = AVT_UNSTRUCTURED_MESH;
	//
	int nblocks = m_number_domain;
	//debug1 << "num of domain "<<m_number_domain<<"\n";
	int block_origin = 0;
	int cell_origin = 0;
	int group_origin = 0;
	int spatial_dimension = 3;
	int topological_dimension = 3;
	double *extents = NULL;
	int * bounds = NULL;
	//add node center 3d mesh
	if (m_external_mesh_providers.begin()->second->provide3DMesh())
	{
		avtMeshMetaData *mmd = new avtMeshMetaData(bounds, extents, mesh_name,
			nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
		a_metaData->Add(mmd);
	}

	//
	// add layered 2d mesh
	mesh_name = m_layer_mesh;
	topological_dimension = 2;

	if (m_external_mesh_providers.begin()->second->provide3DMesh())
	{
		avtMeshMetaData *mmd = new avtMeshMetaData(bounds, extents, mesh_name,
			nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
		a_metaData->Add(mmd);
	}

	//
	// add surface 2d mesh
	//
	mesh_name = m_mesh_2d;

	spatial_dimension = 2;
	topological_dimension = 2;

	avtMeshMetaData *mmd1 = new avtMeshMetaData(bounds, extents, mesh_name,
		nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
	a_metaData->Add(mmd1);

	//add 3d side center point mesh
	mesh_name = m_side_center_point_3d_mesh;

	spatial_dimension = 3;
	topological_dimension = 0;

	if (m_external_mesh_providers.begin()->second->provide3DMesh())
	{
		avtMeshMetaData *mmd = new avtMeshMetaData(bounds, extents, mesh_name,
			nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
		a_metaData->Add(mmd);
	}

	//add 3d  face center point mesh
	mesh_name = m_face_center_point_3d_mesh;

	spatial_dimension = 3;
	topological_dimension = 0;

	if (m_external_mesh_providers.begin()->second->provide3DMesh())
	{
		avtMeshMetaData *mmd = new avtMeshMetaData(bounds, extents, mesh_name,
			nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
		a_metaData->Add(mmd);
	}

	// add 2d side center point mesh
	mesh_name = m_side_center_point_2d_mesh;

	spatial_dimension = 2;
	topological_dimension = 0;

	avtMeshMetaData *mmd2 = new avtMeshMetaData(bounds, extents, mesh_name,
		nblocks, block_origin, cell_origin, group_origin, spatial_dimension, topological_dimension, mt);
	a_metaData->Add(mmd2);

	//
	// add water surface and depth scalar
	string mesh = m_mesh_2d;
	avtScalarMetaData *smd1 = new avtScalarMetaData(m_node_depth_label, mesh, nodeCent);
	a_metaData->Add(smd1);



	// m_var_mesh_map[m_node_surface_label] = mesh;
	m_var_mesh_map[m_node_depth_label] = mesh;

	//add 3D node level label
	string mesh3d = m_mesh_3d;
	avtScalarMetaData *smd2 = new avtScalarMetaData(MeshConstants10::NODE_LEVEL, mesh3d, nodeCent);
	a_metaData->Add(smd2);

	//add 3D element level label
	mesh3d = m_layer_mesh;
	avtScalarMetaData *smd3 = new avtScalarMetaData(MeshConstants10::ELE_LEVEL, mesh3d, zoneCent);
	a_metaData->Add(smd3);

	//add 3D side level label
	mesh3d = m_side_center_point_3d_mesh;
	avtScalarMetaData *smd4 = new avtScalarMetaData(MeshConstants10::SIDE_LEVEL, mesh3d, nodeCent);
	a_metaData->Add(smd4);

	//add 3D layer lable (for prism center data)
	mesh3d = m_mesh_3d;
	avtScalarMetaData *smd5 = new avtScalarMetaData(MeshConstants10::LAYER, mesh3d, zoneCent);
	a_metaData->Add(smd5);


	//add node bottom label
	mesh = m_mesh_2d;
	//a_avtFile->addScalarVarToMetaData(a_metaData,  MeshConstants10::NODE_BOTTOM, mesh, nodeCent);
	//a_avtFile->addScalarVarToMetaData(a_metaData,  MeshConstants10::FACE_BOTTOM, mesh, zoneCent);

	//a_avtFile->addScalarVarToMetaData(a_metaData,  MeshConstants10::EDGE_BOTTOM, mesh, nodeCent);


	PopulateStateMetaData(a_metaData, a_timeState);
	debug1 << "finish populate metadata \n";
#endif
	
}

void    avtMDSCHISMFileFormatImpl::addFaceCenterData(avtDatabaseMetaData * a_metaData,
	SCHISMVar10           * a_varPtr,
	const std::string   & a_varName,
	const std::string   & a_varLabel,
	const avtCentering  & a_center)

{
	// only add face centered  var now 
	string mesh2d = m_mesh_2d;
	string mesh3d = m_mesh_3d;
	std::string level_center = a_varPtr->get_vertical_center();

	if (level_center == MeshConstants10::FULL_LAYER)
	{
		mesh3d = m_layer_mesh;
	}

	avtCentering  faceCent(AVT_ZONECENT);
	// scalar data 2d mesh
	if (a_varPtr->num_dims() <= 2)
	{
		avtScalarMetaData *smd = new avtScalarMetaData(a_varLabel, mesh2d, faceCent);
		a_metaData->Add(smd);
		
		m_var_mesh_map[a_varLabel] = mesh2d;
		m_var_dim[a_varName] = a_varPtr->num_dims();
	}
	// vector data
	else if (SCHISMVarIs3D(a_varPtr))
	{
		if (!(m_external_mesh_providers.begin()->second->provide3DMesh()))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "3D variable " << a_varLabel << " is not supported by a 2D meshprovider\n";
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		if (SCHISMVarIsVector(a_varPtr)) // 3d vector
		{
			// last dim is vector component         
			SCHISMDim10* comDim = a_varPtr->get_dim(3);
			int ncomps = comDim->size();
			int ucomps = (ncomps == 2 ? 3 : ncomps);

			avtVectorMetaData *vmd = new avtVectorMetaData(a_varLabel, mesh3d, faceCent, ucomps);
			a_metaData->Add(vmd);
			m_var_mesh_map[a_varLabel] = mesh3d;
			m_var_dim[a_varName] = 3;

			// also add bottom, surface and depth average state option
			avtVectorMetaData *vmd1 = new avtVectorMetaData(a_varLabel + m_surface_state_suffix, m_mesh_2d, faceCent, ucomps);
			a_metaData->Add(vmd1);
			//debug1<<"add  "<<label+m_surface_state_suffix<<" ";
			avtVectorMetaData *vmd2 = new avtVectorMetaData(a_varLabel + m_bottom_state_suffix, m_mesh_2d, faceCent, ucomps);
			a_metaData->Add(vmd2);
			//debug1<<"add  "<<label+m_bottom_state_suffix<<" ";
			avtVectorMetaData *vmd3 = new avtVectorMetaData(a_varLabel + m_depth_average_suffix, m_mesh_2d, faceCent, ucomps);
			a_metaData->Add(vmd3);
			// debug1<<"add  "<< label+m_depth_average_suffix<<" ";
			m_var_name_label_map[a_varLabel + m_surface_state_suffix] = a_varName;
			m_var_name_label_map[a_varLabel + m_bottom_state_suffix] = a_varName;
			m_var_name_label_map[a_varLabel + m_depth_average_suffix] = a_varName;
			m_var_mesh_map[a_varLabel + m_surface_state_suffix] = m_mesh_2d;
			m_var_mesh_map[a_varLabel + m_bottom_state_suffix] = m_mesh_2d;
			m_var_mesh_map[a_varLabel + m_depth_average_suffix] = m_mesh_2d;
		}
		else //3d scalar
		{
			avtScalarMetaData *smd = new avtScalarMetaData(a_varLabel, mesh3d, faceCent);
			a_metaData->Add(smd);
			m_var_mesh_map[a_varLabel] = mesh3d;

			// also add bottom, surface and depth average state option
			a_metaData->Add( new avtScalarMetaData(a_varLabel + m_surface_state_suffix,m_mesh_2d,faceCent));
			a_metaData->Add(new avtScalarMetaData(a_varLabel + m_bottom_state_suffix, m_mesh_2d, faceCent));
			a_metaData->Add(new avtScalarMetaData(a_varLabel + m_depth_average_suffix, m_mesh_2d, faceCent));
			//debug1<<"add  "<< label+m_depth_average_suffix<<" ";
			// all those surface, botootm and average are based on original data set
			m_var_name_label_map[a_varLabel + m_surface_state_suffix] = a_varName;
			m_var_name_label_map[a_varLabel + m_bottom_state_suffix] = a_varName;
			m_var_name_label_map[a_varLabel + m_depth_average_suffix] = a_varName;

			m_var_mesh_map[a_varLabel + m_surface_state_suffix] = m_mesh_2d;
			m_var_mesh_map[a_varLabel + m_bottom_state_suffix] = m_mesh_2d;
			m_var_mesh_map[a_varLabel + m_depth_average_suffix] = m_mesh_2d;
			m_var_dim[a_varName] = 3;

		}
		//AddScalarVarToMetaData(a_metaData,  MeshConstants10::LEVEL, mesh3d, faceCent);
	}
	else
	{

	}


}
void    avtMDSCHISMFileFormatImpl::addNodeCenterData(avtDatabaseMetaData * a_metaData,
	SCHISMVar10            * a_varPtr,
	const std::string   & a_varName,
	const std::string   & a_varLabel,
	const avtCentering  & a_center)
{

	std::string varName(a_varName);
	std::string label(a_varLabel);
	avtCentering avtCenter(a_center);

	//  scalar var on 2D
	if (a_varPtr->num_dims() <= 2)
	{
		a_metaData->Add(new avtScalarMetaData(label, m_mesh_2d, avtCenter));
		m_var_mesh_map[label] = m_mesh_2d;
		m_var_dim[varName] = a_varPtr->num_dims();
		debug1 << "added 2d scalar:" << label;
	}
	//  vector var on 2D
	else if ((a_varPtr->num_dims() == 3) && (!SCHISMVarIs3D(a_varPtr)))
	{
		SCHISMDim10* comDim = a_varPtr->get_dim(3);
		int ncomps = comDim->size();
		int ucomps = (ncomps == 2 ? 3 : ncomps);
		a_metaData->Add(new avtVectorMetaData(label, m_mesh_2d, avtCenter, ucomps));
		m_var_mesh_map[label] = m_mesh_2d;
		m_var_dim[varName] = 2;

	}
	//  scalar var having layer dim
	else if ((a_varPtr->num_dims() == 3) && (SCHISMVarIs3D(a_varPtr)))
	{
		if (!(m_external_mesh_providers.begin()->second->provide3DMesh()))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "3D variable " << a_varLabel << " is not supported by a 2D meshprovider\n";
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}

		//AddScalarVarToMetaData(a_metaData,  MeshConstants10::LEVEL, m_mesh_3d, avtCenter);
		a_metaData->Add(new avtScalarMetaData(label, m_mesh_3d, avtCenter));
		m_var_mesh_map[label] = m_mesh_3d;
		// also add bottom, surface and depth average state option
		if (!(varName == MeshConstants10::ZCOORD))
		{
			a_metaData->Add(new avtScalarMetaData(
				label + m_surface_state_suffix,
				m_mesh_2d,
				avtCenter));
			debug1 << "add  " << label + m_surface_state_suffix << " ";
			a_metaData->Add(new avtScalarMetaData(
				label + m_bottom_state_suffix,
				m_mesh_2d,
				avtCenter));
			debug1 << "add  " << label + m_bottom_state_suffix << " ";
			a_metaData->Add(new avtScalarMetaData(
				label + m_depth_average_suffix,
				m_mesh_2d,
				avtCenter));
			debug1 << "add  " << label + m_depth_average_suffix << " ";


			// all those surface, botootm and average are based on original data set
			m_var_name_label_map[label + m_surface_state_suffix] = varName;
			m_var_name_label_map[label + m_bottom_state_suffix] = varName;
			m_var_name_label_map[label + m_depth_average_suffix] = varName;

			m_var_mesh_map[label + m_surface_state_suffix] = m_mesh_2d;
			m_var_mesh_map[label + m_bottom_state_suffix] = m_mesh_2d;
			m_var_mesh_map[label + m_depth_average_suffix] = m_mesh_2d;
		}
		m_var_dim[varName] = 3;

	}
	else if ((a_varPtr->num_dims() == 4) && (SCHISMVarIs3D(a_varPtr)))
	{
		if (!(m_external_mesh_providers.begin()->second->provide3DMesh()))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "3D variable " << a_varLabel << " is not supported by a 2D meshprovider\n";
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		// last dim is vector component         
		SCHISMDim10* comDim = a_varPtr->get_dim(3);
		int ncomps = comDim->size();
		int ucomps = (ncomps == 2 ? 3 : ncomps);

		a_metaData->Add(new avtVectorMetaData(label, m_mesh_3d, avtCenter, ucomps));
		//AddScalarVarToMetaData(a_metaData,  MeshConstants10::LEVEL, m_mesh_3d, avtCenter);
		m_var_mesh_map[label] = m_mesh_3d;

		// also add bottom, surface and depth average state option
		a_metaData->Add(new avtVectorMetaData(
			label + m_surface_state_suffix,
			m_mesh_2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_surface_state_suffix << " ";
		a_metaData->Add(new avtVectorMetaData(
			label + m_bottom_state_suffix,
			m_mesh_2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_bottom_state_suffix << " ";
		a_metaData->Add(new avtVectorMetaData(
			label + m_depth_average_suffix,
			m_mesh_2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_depth_average_suffix << " ";
		m_var_name_label_map[label + m_surface_state_suffix] = varName;
		m_var_name_label_map[label + m_bottom_state_suffix] = varName;
		m_var_name_label_map[label + m_depth_average_suffix] = varName;
		m_var_mesh_map[label + m_surface_state_suffix] = m_mesh_2d;
		m_var_mesh_map[label + m_bottom_state_suffix] = m_mesh_2d;
		m_var_mesh_map[label + m_depth_average_suffix] = m_mesh_2d;
		m_var_dim[varName] = 3;

	}
	else
	{

	}
}

void    avtMDSCHISMFileFormatImpl::addSideCenterData(avtDatabaseMetaData * a_metaData,
	SCHISMVar10           * a_varPtr,
	const std::string   & a_varName,
	const std::string   & a_varLabel,
	const avtCentering  & a_center)
{
	// only add centered  var now 
	string mesh2d = m_side_center_point_2d_mesh;
	string mesh3d = m_side_center_point_3d_mesh;
	std::string level_center = a_varPtr->get_vertical_center();

	if (level_center == MeshConstants10::HALF_LAYER)
	{
		mesh3d = m_face_center_point_3d_mesh;
	}

	avtCentering  nodeCent(AVT_NODECENT);
	std::string varName(a_varName);
	std::string label(a_varLabel);
	avtCentering avtCenter(a_center);

	//  scalar var on 2D
	if (a_varPtr->num_dims() <= 2)
	{
		a_metaData->Add(new avtScalarMetaData(label, mesh2d, avtCenter));
		m_var_mesh_map[label] = mesh2d;
		m_var_dim[varName] = a_varPtr->num_dims();
		debug1 << "added 2d scalar:" << label;
	}
	//  vector var on 2D
	else if ((a_varPtr->num_dims() == 3) && (!SCHISMVarIs3D(a_varPtr)))
	{
		SCHISMDim10* comDim = a_varPtr->get_dim(3);
		int ncomps = comDim->size();
		int ucomps = (ncomps == 2 ? 3 : ncomps);
		a_metaData->Add(new avtVectorMetaData(label, mesh2d, avtCenter, ucomps));
		m_var_mesh_map[label] = mesh2d;
		m_var_dim[varName] = 2;

	}
	//  scalar var having layer dim
	else if ((a_varPtr->num_dims() == 3) && (SCHISMVarIs3D(a_varPtr)))
	{

		if (!(m_external_mesh_providers.begin()->second->provide3DMesh()))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "3D variable " << a_varLabel << " is not supported by a 2D meshprovider\n";
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		//AddScalarVarToMetaData(a_metaData,  MeshConstants10::LEVEL, mesh3d, avtCenter);
		a_metaData->Add(new avtScalarMetaData(label, mesh3d, avtCenter));
		m_var_mesh_map[label] = mesh3d;
		// also add bottom, surface and depth average state option
		a_metaData->Add(new avtScalarMetaData(
			label + m_surface_state_suffix,
			mesh2d,
			avtCenter));
		debug1 << "add  " << label + m_surface_state_suffix << " ";
		a_metaData->Add(new avtScalarMetaData(
			label + m_bottom_state_suffix,
			mesh2d,
			avtCenter));
		debug1 << "add  " << label + m_bottom_state_suffix << " ";
		a_metaData->Add(new avtScalarMetaData(
			label + m_depth_average_suffix,
			mesh2d,
			avtCenter));
		debug1 << "add  " << label + m_depth_average_suffix << " ";
		// all those surface, botootm and average are based on original data set
		m_var_name_label_map[label + m_surface_state_suffix] = varName;
		m_var_name_label_map[label + m_bottom_state_suffix] = varName;
		m_var_name_label_map[label + m_depth_average_suffix] = varName;

		m_var_mesh_map[label + m_surface_state_suffix] = mesh2d;
		m_var_mesh_map[label + m_bottom_state_suffix] = mesh2d;
		m_var_mesh_map[label + m_depth_average_suffix] = mesh2d;
		m_var_dim[varName] = 3;

	}
	else if ((a_varPtr->num_dims() == 4) && (SCHISMVarIs3D(a_varPtr)))
	{
		if (!(m_external_mesh_providers.begin()->second->provide3DMesh()))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "3D variable " << a_varLabel << " is not supported by a 2D meshprovider\n";
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		// last dim is vector component         
		SCHISMDim10* comDim = a_varPtr->get_dim(3);
		int ncomps = comDim->size();
		int ucomps = (ncomps == 2 ? 3 : ncomps);

		a_metaData->Add(new avtVectorMetaData(label, mesh3d, avtCenter, ucomps));
		//AddScalarVarToMetaData(a_metaData,  MeshConstants10::LEVEL, mesh3d, avtCenter);
		m_var_mesh_map[label] = mesh3d;

		// also add bottom, surface and depth average state option
		a_metaData->Add(new avtVectorMetaData(
			label + m_surface_state_suffix,
			mesh2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_surface_state_suffix << " ";
		a_metaData->Add(new avtVectorMetaData(
			label + m_bottom_state_suffix,
			mesh2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_bottom_state_suffix << " ";
		a_metaData->Add(new avtVectorMetaData(
			label + m_depth_average_suffix,
			mesh2d,
			avtCenter,
			ucomps));
		debug1 << "add  " << label + m_depth_average_suffix << " ";
		m_var_name_label_map[label + m_surface_state_suffix] = varName;
		m_var_name_label_map[label + m_bottom_state_suffix] = varName;
		m_var_name_label_map[label + m_depth_average_suffix] = varName;
		m_var_mesh_map[label + m_surface_state_suffix] = mesh2d;
		m_var_mesh_map[label + m_bottom_state_suffix] = mesh2d;
		m_var_mesh_map[label + m_depth_average_suffix] = mesh2d;
		m_var_dim[varName] = 3;
	}
	else
	{
	}

}

MDSchismOutput * avtMDSCHISMFileFormatImpl::get_a_data_file()
{
	return m_data_files.begin()->second;
}
MDSCHISMMeshProvider * avtMDSCHISMFileFormatImpl::get_a_mesh_provider()
{
	return m_external_mesh_providers.begin()->second;
}

// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::PopulateStateMetaData
//
//  Purpose:
//      Scan the data file and find out all the state variables, 
//      store the name, mesh, and centering information into
//      the input metadata set.
//
//  Programmer: qshu
//  Creation:   Wed Aug 30 08:11:04 PDT 2012
//
// ****************************************************************************

void avtMDSCHISMFileFormatImpl::PopulateStateMetaData(avtDatabaseMetaData * a_metaData,
	int                   a_timeState)
{
	int numVar = m_data_files.begin()->second->num_vars();
	//debug1 << "get vars " << numVar << endl;
	for (int iVar = 0; iVar < numVar; iVar++)
	{
		//debug1 << iVar;
		SCHISMVar10*  varPtr = m_data_files.begin()->second->get_var(iVar);
		std::string varName = varPtr->name();
		
		if (m_data_files.begin()->second->none_data_var(varName))
		{
			debug1 << varName << " is skipped\n";
			continue;
		}
		debug1 << iVar << " " << varPtr->num_dims() << " " << varName << endl;
		if (!(varPtr->is_defined_over_grid()))
		{
			debug1 << varName << " passed \n";
			continue;
		}

		if (varPtr->is_SCHISM_mesh_parameter())
		{
			debug1 << varName << "passed 2\n";
			continue;
		}

		//debug1 << varName<<"to be added" << endl;

		std::string  location(NODE);
		avtCentering avtCenter(AVT_NODECENT);


		if ((varName == m_node_surface) || (varName == m_node_depth))
		{
			continue;
		}
		std::string  label;
		label = varName;
		// this dic make it easy to find out data set for a visit plot variable
		m_var_name_label_map[label] = varName;

		// handle different for face and node center data
		location = varPtr->get_horizontal_center();

		std::string vertical_center = varPtr->get_vertical_center();

		m_var_horizontal_center_map[varName] = location;
		m_var_vertical_center_map[varName] = vertical_center;
		
		if (location == FACE)
		{
			//debug1 << " begin add face var meta\n";
			addFaceCenterData(a_metaData, varPtr, varName, label, avtCenter);
		}
		else if (location == NODE)
		{
			//debug1 << " begin add node var meta\n";
			addNodeCenterData(a_metaData, varPtr, varName, label, avtCenter);
		}
		else if (location == SIDE)
		{
			//debug1 << " begin add side var meta\n";
			addSideCenterData(a_metaData, varPtr, varName, label, avtCenter);
		}
		// omit unkown center data
		else
		{
			continue;
		}
	}
}


void   avtMDSCHISMFileFormatImpl::create2DUnstructuredMesh(vtkUnstructuredGrid *a_uGrid,
	long                 *a_meshEle,
	const  int          &a_domainID,
	const  int          &a_timeState)
{
	long   numNodes = m_num_mesh_nodes[a_domainID];
	vtkPoints *points = vtkPoints::New();
	points->SetDataTypeToDouble();
	points->GetData()->SetNumberOfComponents(3);
	points->SetNumberOfPoints(numNodes);

	points->SetNumberOfPoints(numNodes);
	double * pointPtr = (double *)points->GetVoidPointer(0);

	if (!m_external_mesh_providers[a_domainID]->fillPointCoord2D(pointPtr, a_timeState))
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve faces nodes coord at step " << a_timeState;
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}
	/*vtkUnsignedCharArray *ghost_zones = vtkUnsignedCharArray::New();
	ghost_zones->SetName("avtGhostZones");
	long ncells = m_num_mesh_faces[a_domainID];
	ghost_zones->SetNumberOfTuples(ncells);
	unsigned char *gzp = ghost_zones->GetPointer(0);
	for (long i = 0; i < ncells; i++)
		gzp[i] = 0;
	unsigned char val = 1;
	avtGhostData::AddGhostZoneType(val, DUPLICATED_ZONE_INTERNAL_TO_PROBLEM);
	long ncell_not_ghost = m_external_mesh_providers[a_domainID]->m_number_element_no_ghost;
	for (long i = ncell_not_ghost; i < ncells; i++)
		gzp[i] = val;*/
	unsigned char val2 = 1;
	avtGhostData::AddGhostNodeType(val2, DUPLICATED_NODE);
	vtkUnsignedCharArray *ghost_nodes = vtkUnsignedCharArray::New();
	ghost_nodes->SetName("avtGhostNodes");
	long npts = m_num_mesh_nodes[a_domainID];
	ghost_nodes->SetNumberOfTuples(npts);
	unsigned char *gnp = ghost_nodes->GetPointer(0);
	for (long i = 0; i < npts; i++)
	{
		long global_id = m_external_mesh_providers[a_domainID]->m_local_node_id_to_global_id[i];
		int num_domain = m_node_num_domain[global_id - 1];
		if (num_domain>1)
			gnp[i] = val2;
		else
		    gnp[i] = 0;
	}
	
	
	a_uGrid->SetPoints(points);
	points->Delete();
	a_uGrid->Allocate(m_num_mesh_faces[a_domainID]);
	
	long *  nodePtrTemp = a_meshEle;
	load_ele_dry_wet(a_timeState,a_domainID);
	for (long iCell = 0; iCell < m_num_mesh_faces[a_domainID]; ++iCell)
	{
		int numberOfNodeInCell = *nodePtrTemp;

		//if (!(m_ele_dry_wet[iCell]))
		//{
		if (numberOfNodeInCell == 3)
		{
			vtkIdType verts[3];
			for (int iNode = 0; iNode < 3; ++iNode)
			{
				verts[iNode] = nodePtrTemp[iNode + 1] - 1;

			}
			nodePtrTemp += (MeshConstants10::MAX_NUM_NODE_PER_CELL + 1);

			a_uGrid->InsertNextCell(VTK_TRIANGLE, 3, verts);
		}
		else if (numberOfNodeInCell == 4)
		{
			vtkIdType verts[4];
			for (int iNode = 0; iNode < 4; ++iNode)
			{
				verts[iNode] = nodePtrTemp[iNode + 1] - 1;

			}
			nodePtrTemp += (MeshConstants10::MAX_NUM_NODE_PER_CELL + 1);

			a_uGrid->InsertNextCell(VTK_QUAD, 4, verts);
		}
		else
		{
			stringstream msgStream(stringstream::out);
			msgStream << "invalid cell type with number of nodes: " << numberOfNodeInCell;
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		//}

	}
	

	//a_uGrid->GetCellData()->AddArray(ghost_zones);
	a_uGrid->GetCellData()->AddArray(ghost_nodes);
	
}


void   avtMDSCHISMFileFormatImpl::create2DUnstructuredMeshNoDryWet(vtkUnstructuredGrid *a_uGrid,
	long                 *a_meshEle,
	const  int          &a_domainID,
	const  int          &a_timeState)
{
	long   numNodes = m_num_mesh_nodes[a_domainID];
	vtkPoints *points = vtkPoints::New();
	points->SetDataTypeToDouble();
	points->GetData()->SetNumberOfComponents(3);
	points->SetNumberOfPoints(numNodes);
	double * pointPtr = (double *)points->GetVoidPointer(0);

	if (!m_external_mesh_providers[a_domainID]->fillPointCoord2D(pointPtr, a_timeState))
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve faces nodes coord at step " << a_timeState;
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}


	a_uGrid->SetPoints(points);
	points->Delete();
	a_uGrid->Allocate(m_num_mesh_faces[a_domainID]);

	long *  nodePtrTemp = a_meshEle;

	for (long iCell = 0; iCell < m_num_mesh_faces[a_domainID]; ++iCell)
	{
		int numberOfNodeInCell = *nodePtrTemp;


		if (numberOfNodeInCell == 3)
		{
			vtkIdType verts[3];
			for (int iNode = 0; iNode < 3; ++iNode)
			{
				verts[iNode] = nodePtrTemp[iNode + 1] - 1;

			}
			nodePtrTemp += (MeshConstants10::MAX_NUM_NODE_PER_CELL + 1);

			a_uGrid->InsertNextCell(VTK_TRIANGLE, 3, verts);
		}
		else if (numberOfNodeInCell == 4)
		{
			vtkIdType verts[4];
			for (int iNode = 0; iNode < 4; ++iNode)
			{
				verts[iNode] = nodePtrTemp[iNode + 1] - 1;

			}
			nodePtrTemp += (MeshConstants10::MAX_NUM_NODE_PER_CELL + 1);

			a_uGrid->InsertNextCell(VTK_QUAD, 4, verts);
		}
		else
		{
			stringstream msgStream(stringstream::out);
			msgStream << "invalid cell type with number of nodes: " << numberOfNodeInCell;
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}


	}

}


void   avtMDSCHISMFileFormatImpl::createLayerMesh(vtkUnstructuredGrid *a_uGrid,
	long                 *a_meshEle,
	long                *a_2DPointto3DPoints,
	const  int          &a_domainID,
	const  int          &a_timeState)
{
	vtkPoints *points = vtkPoints::New();
	points->SetDataTypeToDouble();
	points->GetData()->SetNumberOfComponents(3);
	points->SetNumberOfPoints(m_total_valid_3D_point[a_domainID]);
	double * pointPtr = (double *)points->GetVoidPointer(0);
	//debug only

	if (!m_external_mesh_providers[a_domainID]->fillPointCoord3D(pointPtr, a_timeState))
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve faces nodes coord at step " << a_timeState;
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}

	a_uGrid->SetPoints(points);
	points->Delete();
	a_uGrid->Allocate(m_num_mesh_faces[a_domainID]*m_num_layers);


	int * kbe = m_kbp_ele[a_domainID];
	debug1 << "test layer mesh ele\n";
	for (int iLayer = 0; iLayer < m_num_layers; iLayer++)
	{

		for (int iCell = 0; iCell < m_num_mesh_faces[a_domainID]; ++iCell)
		{

			if (iLayer >= (std::max(1, kbe[iCell]) - 1))
			{

				int numberOfNodeInCell = a_meshEle[iCell*(NODESPERELE + 1)];

				if (numberOfNodeInCell == 3)
				{
					vtkIdType verts[3];

					for (int i = 0; i < 3; i++)
					{
						int p = a_meshEle[iCell*(NODESPERELE + 1) + i + 1] - 1;
						int p3d = a_2DPointto3DPoints[p*m_num_layers + iLayer];
						int valid_bottom = std::max(1, m_kbp_node[a_domainID][p]) - 1;
						if (iLayer < valid_bottom)
						{
							p3d = a_2DPointto3DPoints[p*m_num_layers + valid_bottom];
						}

						verts[i] = p3d;
					}

					a_uGrid->InsertNextCell(VTK_TRIANGLE, 3, verts);
				}
				else if (numberOfNodeInCell == 4)
				{
					vtkIdType verts[4];
					for (int i = 0; i < 4; i++)
					{
						int p = a_meshEle[iCell*(NODESPERELE + 1) + i + 1] - 1;
						int p3d = a_2DPointto3DPoints[p*m_num_layers + iLayer];
						int valid_bottom = std::max(1, m_kbp_node[a_domainID][p]) - 1;
						if (iLayer < valid_bottom)
						{
							p3d = a_2DPointto3DPoints[p*m_num_layers + valid_bottom];
						}

						verts[i] = p3d;
					}
					a_uGrid->InsertNextCell(VTK_QUAD, 4, verts);
				}
				else
				{
					stringstream msgStream(stringstream::out);
					msgStream << "invalid cell type with number of nodes: " << numberOfNodeInCell;
					EXCEPTION1(InvalidVariableException, msgStream.str());
				}
			}

		}

	}


}


void   avtMDSCHISMFileFormatImpl::create3DUnstructuredMesh(vtkUnstructuredGrid *a_uGrid,
	long                 *a_meshEle,
	long                 *a_2DPointto3DPoints,
	const  int          &a_domainID,
	const  int          &a_timeState)
{
	vtkPoints *points = vtkPoints::New();
	points->SetDataTypeToDouble();
	points->GetData()->SetNumberOfComponents(3);
	points->SetNumberOfPoints(m_total_valid_3D_point[a_domainID]);
	//debug1 << "total valid 3d pts: " << m_total_valid_3D_point[a_domainID] << "\n";
	double * pointPtr = (double *)points->GetVoidPointer(0);

	if (!m_external_mesh_providers[a_domainID]->fillPointCoord3D(pointPtr, a_timeState))
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve faces nodes coord at step " << a_timeState;
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}

	unsigned char val2 = 1;
	avtGhostData::AddGhostNodeType(val2, DUPLICATED_NODE);
	vtkUnsignedCharArray *ghost_nodes = vtkUnsignedCharArray::New();
	ghost_nodes->SetName("avtGhostNodes");
	long npts = m_total_valid_3D_point[a_domainID];
	ghost_nodes->SetNumberOfTuples(npts);
	unsigned char *gnp = ghost_nodes->GetPointer(0);
	for (long i = 0; i < npts; i++)
	{
		long ii = m_3d_node_to_2d_node[a_domainID][i];
		long global_id = m_external_mesh_providers[a_domainID]->m_local_node_id_to_global_id[ii];
		int num_domain = m_node_num_domain[global_id - 1];
		if (num_domain > 1)
			gnp[i] = val2;
		else
			gnp[i] = 0;
	}

	debug1 << "finish compute and cahce z\n";
	a_uGrid->SetPoints(points);
	points->Delete();
	a_uGrid->Allocate(m_num_mesh_faces[a_domainID]*(m_num_layers - 1));
	debug1 << "total 2d cell num " << m_num_mesh_faces[a_domainID] << "\n";
	long *  nodePtrTemp = a_meshEle;


	m_tri_wedge = 0;
	m_tri_pyramid = 0;
	m_tri_tetra = 0;

	m_quad_hexhedron = 0;
	m_quad_wedge = 0;
	m_quad_pyramid = 0;
	load_ele_dry_wet(a_timeState,a_domainID);

	bool *prism_bottom_set = new bool[m_num_mesh_faces[a_domainID]];
	for (int iCell = 0; iCell < m_num_mesh_faces[a_domainID]; ++iCell)
	{
		prism_bottom_set[iCell] = false;
	}


	for (int iLayer = 0; iLayer < m_num_layers - 1; iLayer++)
	{
		nodePtrTemp = a_meshEle;
		for (int iCell = 0; iCell < m_num_mesh_faces[a_domainID]; ++iCell)
		{
			{
				nodePtrTemp = a_meshEle + (MeshConstants10::MAX_NUM_NODE_PER_CELL + 1)*iCell;
				int numberOfNodeInCell = *nodePtrTemp;

				long validTopNode[MeshConstants10::MAX_NUM_NODE_PER_CELL];
				long validBottomNode[MeshConstants10::MAX_NUM_NODE_PER_CELL];
				int validTopNodeNum = 0;
				int validBottomNodeNum = 0;

				validTopBottomNode(validTopNodeNum,
					validBottomNodeNum,
					validTopNode,
					validBottomNode,
					iLayer,
					nodePtrTemp,
					m_kbp_node[a_domainID]);

				if (!(prism_bottom_set[iCell]))
				{
					if (validBottomNodeNum > 0)
					{
						m_kbp_prism[a_domainID][iCell] = iLayer+1; //level starts from 1 schism tradition
						prism_bottom_set[iCell] = true;
					}
				}


				if (numberOfNodeInCell == 3)
				{
					insertTriangle3DCell(a_uGrid,
						validTopNodeNum,
						validBottomNodeNum,
						validTopNode,
						validBottomNode,
						nodePtrTemp,
						a_2DPointto3DPoints,
						m_kbp_node[a_domainID],
						iCell,
						iLayer,
						m_num_layers);
					//move pointer to next element
					//nodePtrTemp += (MeshConstants10::MAX_NUM_NODE_PER_CELL+1);
				}
				else if (numberOfNodeInCell == 4)
				{
					insertQuad3DCell(a_uGrid,
						validTopNodeNum,
						validBottomNodeNum,
						validTopNode,
						validBottomNode,
						nodePtrTemp,
						a_2DPointto3DPoints,
						m_kbp_node[a_domainID],
						iCell,
						iLayer,
						m_num_layers);
					// nodePtrTemp += (MeshConstants10::MAX_NUM_NODE_PER_CELL+1);
				}

				else
				{
					//omit
				}

			}
		}

	}
	m_data_files[a_domainID]->set_prism_bottom(a_timeState, m_kbp_prism[a_domainID]);
	a_uGrid->GetCellData()->AddArray(ghost_nodes);
	delete [] prism_bottom_set;
	debug1 << " tri_wedge " << m_tri_wedge << " tri_pyramid " << m_tri_pyramid << " tri_tetra " << m_tri_tetra;

	debug1 << " quad_hexhedron " << m_quad_hexhedron << " quad_wedge " << m_quad_wedge << " quad_pyramid " << m_quad_pyramid << "\n";
}

void    avtMDSCHISMFileFormatImpl::create2DPointMesh(vtkUnstructuredGrid *a_uGrid,
	long                 *a_meshEle,
	const  int          &a_domainID,
	const  int          &a_timeState)
{
	long   numNodes = m_num_mesh_edges[a_domainID];
	vtkPoints *points = vtkPoints::New();
	points->SetDataTypeToDouble();
	points->GetData()->SetNumberOfComponents(3);
	points->SetNumberOfPoints(numNodes);
	double * pointPtr = (double *)points->GetVoidPointer(0);

	if (!m_external_mesh_providers[a_domainID]->fillSideCenterCoord2D(pointPtr, a_timeState))
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve edge center coord at step " << a_timeState;
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}


	a_uGrid->SetPoints(points);
	points->Delete();
	a_uGrid->Allocate(numNodes);
	vtkIdType onevertex;
	for (long i = 0; i < numNodes; ++i)
	{
		onevertex = i;
		a_uGrid->InsertNextCell(VTK_VERTEX, 1, &onevertex);
	}

}

void   avtMDSCHISMFileFormatImpl::create3DPointMesh(vtkUnstructuredGrid *a_uGrid,
	long            *a_meshEle,
	const  int     &a_domainID,
	const  int     &a_timeState)
{
	long   numNodes = m_total_valid_3D_side[a_domainID];
	vtkPoints *points = vtkPoints::New();
	points->SetDataTypeToDouble();
	points->GetData()->SetNumberOfComponents(3);
	points->SetNumberOfPoints(numNodes);
	double * pointPtr = (double *)points->GetVoidPointer(0);

	if (!m_external_mesh_providers[a_domainID]->fillSideCenterCoord3D(pointPtr, a_timeState))
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve edge center coord at step " << a_timeState;
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}

	a_uGrid->SetPoints(points);
	points->Delete();
	a_uGrid->Allocate(numNodes);
	vtkIdType onevertex;
	for (long i = 0; i < numNodes; ++i)
	{
		onevertex = i;
		a_uGrid->InsertNextCell(VTK_VERTEX, 1, &onevertex);
	}

}

// this is the mesh consists of center at 3d prism side face, used for var like flux 
void   avtMDSCHISMFileFormatImpl::create3DPointFaceMesh(vtkUnstructuredGrid *a_uGrid,
	long               *a_meshEle,
	const  int         &a_domainID,
	const  int         &a_timeState)
{
	int   numNodes = m_total_valid_3D_side[a_domainID] - (m_external_mesh_providers[a_domainID]->numberOfSide());
	vtkPoints *points = vtkPoints::New();
	points->SetDataTypeToDouble();
	points->GetData()->SetNumberOfComponents(3);
	points->SetNumberOfPoints(numNodes);
	double * pointPtr = (double *)points->GetVoidPointer(0);

	if (!m_external_mesh_providers[a_domainID]->fillSideFaceCenterCoord3D(pointPtr, a_timeState))
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve edge center coord at step " << a_timeState;
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}

	a_uGrid->SetPoints(points);
	points->Delete();
	a_uGrid->Allocate(numNodes);
	vtkIdType onevertex;
	for (int i = 0; i < numNodes; ++i)
	{
		onevertex = i;
		a_uGrid->InsertNextCell(VTK_VERTEX, 1, &onevertex);
	}

}

// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:+
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      mesh_name    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Wed Mar 13 09:13:49 PDT 2013
//
// ****************************************************************************

vtkDataSet *
avtMDSCHISMFileFormatImpl::GetMesh(int a_timeState, int a_domainID, avtMDSCHISMFileFormat * a_avtFile, const char *mesh_name)
{
	debug1 << "begin get mesh\n";
	int   nDims = 3;
	long   numNodes = m_num_mesh_nodes[a_domainID];
	long   numCells = m_num_mesh_faces[a_domainID];

	time_t startTicks = clock();

	int   domainID = a_domainID;
	int   timeState = a_timeState;
	std::string material("all");
	std::string cacheMeshID(mesh_name);
	cacheMeshID += m_data_file;
	cacheMeshID += std::to_string(domainID);
	debug1 << " try to find " << cacheMeshID << " in cache\n";
	//vtkObject * cachedMesh = NULL;
	//cachedMesh = (a_avtFile->get_cache())->GetVTKObject(cacheMeshID.c_str(),
	//	avtVariableCache::DATASET_NAME,
	//	a_timeState,
	//	domainID,
	//	material.c_str());



	//if ((cachedMesh != NULL) && (m_mesh_is_static))
	//{
	//	vtkUnstructuredGrid *uGrid = (vtkUnstructuredGrid *)cachedMesh;
	//	uGrid->Register(NULL);
	//	debug1 << "get " << mesh_name << " from cache" << endl;
	//	time_t endTicks = clock();
	//	debug1 << "time used in building mesh :" << endTicks - startTicks << endl;
	//	updateMeshZCoordinates(uGrid->GetPoints(),
	//		a_timeState,
	//		a_domainID,
	//		mesh_name);
	//	return uGrid;
	//}
	debug1 << mesh_name << " not in cache/or not static. Load from data." << endl;

	// get face nodes
	int    numNodesPerFace = NODESPERELE;
	long *  faceNodesPtr = new long[m_num_mesh_faces[a_domainID]*(numNodesPerFace + 1)];
	
	//m_external_mesh_providers[a_domainID]->fillMeshElement(faceNodesPtr);
	if (!m_external_mesh_providers[a_domainID]->fillMeshElement(faceNodesPtr))
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve domain faces nodes at step " << a_timeState;
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}

	//debug1 << "get face nodes"<<faceNodesPtr[500]<<" "<<faceNodesPtr[501]<<" "<<faceNodesPtr[502]<<"\n";

	if (!m_mesh_is_static)
	{
		if (m_external_mesh_providers[a_domainID]->update_bottom_layer(a_timeState))
		{
			this->load_bottom(a_timeState,a_domainID);
		}
		//else if (cachedMesh != NULL)
		//{
		//	vtkUnstructuredGrid *uGrid = (vtkUnstructuredGrid *)cachedMesh;
		//	uGrid->Register(NULL);
		//	updateMeshZCoordinates(uGrid->GetPoints(),
		//		a_timeState,
		//		a_domainID,
		//		mesh_name);
		//	return uGrid;
		//}
	}


	long *  m2DPointto3DPoints = new long[m_num_mesh_nodes[a_domainID]*m_num_layers];

	for (long i = 0; i < m_num_mesh_nodes[a_domainID]*m_num_layers; i++)
	{
		m2DPointto3DPoints[i] = MeshConstants10::INVALID_NUM;
	}

	int Index = 0;
	for (int iLayer = 0; iLayer < m_num_layers; iLayer++)
	{
		for (long iNode = 0; iNode < m_num_mesh_nodes[a_domainID]; iNode++)
		{
			int bottomLayer = m_kbp_node[a_domainID][iNode];
			if (bottomLayer <= (iLayer + 1))
			{
				m2DPointto3DPoints[iLayer + iNode * m_num_layers] = Index;
				Index++;
			}

		}

	}

	debug1 << "built mesh " << mesh_name << endl;
	vtkUnstructuredGrid *uGrid = vtkUnstructuredGrid::New();
	if (!strcmp(mesh_name, m_mesh_2d.c_str()))
	{
		create2DUnstructuredMesh(uGrid, faceNodesPtr, a_domainID,a_timeState);
	}
	else if (!strcmp(mesh_name, m_mesh_2d_no_wet_dry.c_str()))
	{
		create2DUnstructuredMeshNoDryWet(uGrid, faceNodesPtr, a_domainID,a_timeState);
	}
	else if (!strcmp(mesh_name, m_layer_mesh.c_str()))
	{
		createLayerMesh(uGrid, faceNodesPtr, m2DPointto3DPoints, a_domainID,a_timeState);
	}
	else if (!strcmp(mesh_name, m_mesh_3d.c_str()))
	{
		create3DUnstructuredMesh(uGrid, faceNodesPtr, m2DPointto3DPoints,a_domainID, a_timeState);
	}
	else if (!strcmp(mesh_name, m_side_center_point_2d_mesh.c_str()))
	{
		create2DPointMesh(uGrid, faceNodesPtr, a_domainID,a_timeState);
	}
	else if (!strcmp(mesh_name, m_side_center_point_3d_mesh.c_str()))
	{

		create3DPointMesh(uGrid, faceNodesPtr, a_domainID,a_timeState);
	}
	else if (!strcmp(mesh_name, m_face_center_point_3d_mesh.c_str()))
	{

		create3DPointFaceMesh(uGrid, faceNodesPtr, domainID,a_timeState);
	}
	else
	{
		EXCEPTION1(InvalidVariableException, mesh_name);
	}

	/*(a_avtFile->get_cache())->CacheVTKObject(cacheMeshID.c_str(),
		avtVariableCache::DATASET_NAME,
		a_timeState,
		domainID,
		material.c_str(),
		uGrid);*/

	time_t endTicks = clock();
	debug1 << "time used in building mesh :" << endTicks - startTicks << endl;
	delete    faceNodesPtr;
	delete    m2DPointto3DPoints;
	debug1 << "finish building mesh\n";
	return    uGrid;
}









void   avtMDSCHISMFileFormatImpl::insertTriangle3DCell(vtkUnstructuredGrid * a_uGrid,
	const int           & a_validTopNodeNum,
	const int           & a_validBottomNodeNum,
	long                 * a_validTopNode,
	long                 * a_validBottomNode,
	long                 * a_faceNodePtr,
	long                 * a_2DPointto3DPoints,
	int                   *     a_kbp_node,
	const long           & a_Cell,
	const int           & a_layerID,
	const int           & a_num_layers)
{
	if ((a_validBottomNodeNum == 0) || ((a_validBottomNodeNum == 2) && (a_validTopNodeNum == 2))) // no 3d cell at all
	{
		return;
	}


	if (a_validBottomNodeNum == 3) // this is a wedge 
	{
		vtkIdType verts[2 * 3];
		//first add bottom face node
		for (int iNode = 0; iNode < 3; ++iNode)
		{
			long p = a_validBottomNode[iNode];
			verts[iNode] = a_2DPointto3DPoints[p*a_num_layers + a_layerID];
		}
		//then add top face node
		for (int iNode = 3; iNode < 2 * 3; ++iNode)
		{
			long p = a_validTopNode[iNode - 3];
			verts[iNode] = a_2DPointto3DPoints[p*a_num_layers + a_layerID + 1];
		}
		a_uGrid->InsertNextCell(VTK_WEDGE, 2 * 3, verts);
		m_tri_wedge++;
	}
	else if (a_validBottomNodeNum == 2)// bottom have two node, this is a pyramid
	{
		vtkIdType verts[5];

		long p1 = a_validBottomNode[0];
		long p2 = a_validBottomNode[1];

		verts[0] = a_2DPointto3DPoints[p1*a_num_layers + a_layerID];
		verts[1] = a_2DPointto3DPoints[p2*a_num_layers + a_layerID];
		verts[2] = a_2DPointto3DPoints[p2*a_num_layers + a_layerID + 1];
		verts[3] = a_2DPointto3DPoints[p1*a_num_layers + a_layerID + 1];
		long p3 = MeshConstants10::INVALID_NUM;

		for (int i = 0; i < 3; i++)
		{
			if (((a_faceNodePtr[i + 1] - 1) != p1) && ((a_faceNodePtr[i + 1] - 1) != p2))
			{
				p3 = a_faceNodePtr[i + 1] - 1;
			}
		}


		if (p3 == MeshConstants10::INVALID_NUM)
		{
			stringstream msgStream(stringstream::out);
			msgStream << "fail to get pyramid apex for cell " << a_Cell << "on layer " << a_layerID;
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		int valid_bottom = std::max(1, a_kbp_node[p3]) - 1;
		if ((a_layerID + 1) < valid_bottom)
		{
			verts[4] = a_2DPointto3DPoints[p3*a_num_layers + valid_bottom];
		}
		else
		{
			verts[4] = a_2DPointto3DPoints[p3*a_num_layers + a_layerID + 1];
		}
		a_uGrid->InsertNextCell(VTK_PYRAMID, 5, verts);
		m_tri_pyramid++;
	}
	else //tetra
	{
		vtkIdType verts[4];
		long p4 = a_validBottomNode[0];


		for (int i = 0; i < 3; i++)
		{

			long p = a_faceNodePtr[i + 1] - 1;
			int valid_bottom = std::max(1, a_kbp_node[p]) - 1;
			if ((a_layerID + 1) < valid_bottom)
			{
				verts[i] = a_2DPointto3DPoints[p*a_num_layers + valid_bottom];
			}
			else
			{
				verts[i] = a_2DPointto3DPoints[p*a_num_layers + a_layerID + 1];
			}
		}

		verts[3] = a_2DPointto3DPoints[p4*a_num_layers + a_layerID];
		a_uGrid->InsertNextCell(VTK_TETRA, 4, verts);
		m_tri_tetra++;
	}

}

bool     avtMDSCHISMFileFormatImpl::fourPointsCoplanar(double p1[3],
	double p2[3],
	double p3[3],
	double p4[3])
{
	double v1[3], v2[3], v3[3];

	for (int i = 0; i < 3; i++)
	{
		v1[i] = p2[i] - p1[i];
		v2[i] = p3[i] - p1[i];
		v3[i] = p4[i] - p1[i];
	}

	//double normal[3];
	//vtkMath::Cross(v1,v2,normal);
	//double dotvalue = vtkMath::Dot(normal,v3);

	double dotvalue = v1[0] * (v2[1] * v3[2] - v2[2] * v3[1]) - v1[1] * (v2[0] * v3[2] - v2[2] * v3[0]) + v1[2] * (v2[0] * v3[1] - v2[1] * v3[0]);

	if (abs(dotvalue) < 1e-6)
	{
		return true;
	}
	else
	{
		return false;
	}

}


void    avtMDSCHISMFileFormatImpl::insertQuad3DCell(vtkUnstructuredGrid *  a_uGrid,
	const int           & a_validTopNodeNum,
	const int           & a_validBottomNodeNum,
	long                * a_validTopNode,
	long                 * a_validBottomNode,
	long                 * a_faceNodePtr,
	long                 * a_2DPointto3DPoints,
	int                   *     a_kbp_node,
	const long           & a_Cell,
	const int           & a_layerID,
	const int           & a_num_layers)

{


	if (a_validBottomNodeNum == 0)
	{
		return;
	}


	vtkIdType p_t[4];
	//debug1<<"valid 2d top ";
	for (int inode = 0; inode < 4; inode++)
	{
		long v_t = a_validTopNode[inode];

		// layer in m_kbp_node starts from 1, a_layerID starts from 0
		int bottomLayer = a_kbp_node[v_t];
		//debug1<<v_t<<" bottom "<<bottomLayer<<" ";
		if ((a_layerID + 2) >= bottomLayer)
		{
			p_t[inode] = a_2DPointto3DPoints[v_t*a_num_layers + a_layerID + 1];
		}
		else
		{
			p_t[inode] = a_2DPointto3DPoints[v_t*a_num_layers + bottomLayer - 1];
		}
	}



	//debug1<< "valid top 3d node "<<p_t[0]<<" "<<p_t[1]<<" "<<p_t[2]<<" "<<p_t[3]<<"\n";

	double coord1[3], coord2[3], coord3[3], coord4[3];

	double* coord_ptr = a_uGrid->GetPoint(p_t[0]);
	coord1[0] = coord_ptr[0];
	coord1[1] = coord_ptr[1];
	coord1[2] = coord_ptr[2];

	coord_ptr = a_uGrid->GetPoint(p_t[1]);
	coord2[0] = coord_ptr[0];
	coord2[1] = coord_ptr[1];
	coord2[2] = coord_ptr[2];

	coord_ptr = a_uGrid->GetPoint(p_t[2]);
	coord3[0] = coord_ptr[0];
	coord3[1] = coord_ptr[1];
	coord3[2] = coord_ptr[2];

	coord_ptr = a_uGrid->GetPoint(p_t[3]);
	coord4[0] = coord_ptr[0];
	coord4[1] = coord_ptr[1];
	coord4[2] = coord_ptr[2];

	bool topCoplane = fourPointsCoplanar(coord1,
		coord2,
		coord3,
		coord4);

	//debug1<<"top Coplane "<<topCoplane<<" "<<a_validBottomNodeNum<<"\n";

	if ((a_validBottomNodeNum == 4))
	{

		// find out if bottom four points are coplaner
		long v1 = a_validBottomNode[0];
		long v2 = a_validBottomNode[1];
		long v3 = a_validBottomNode[2];
		long v4 = a_validBottomNode[3];

		vtkIdType p1 = a_2DPointto3DPoints[v1*a_num_layers + a_layerID];
		vtkIdType p2 = a_2DPointto3DPoints[v2*a_num_layers + a_layerID];
		vtkIdType p3 = a_2DPointto3DPoints[v3*a_num_layers + a_layerID];
		vtkIdType p4 = a_2DPointto3DPoints[v4*a_num_layers + a_layerID];

		//debug1<<"valid bottom 3d node "<<p1<<" "<<p2<<" "<<p3<<" "<<p4<<"\n";


		coord_ptr = a_uGrid->GetPoint(p1);
		coord1[0] = coord_ptr[0];
		coord1[1] = coord_ptr[1];
		coord1[2] = coord_ptr[2];

		coord_ptr = a_uGrid->GetPoint(p2);
		coord2[0] = coord_ptr[0];
		coord2[1] = coord_ptr[1];
		coord2[2] = coord_ptr[2];

		coord_ptr = a_uGrid->GetPoint(p3);
		coord3[0] = coord_ptr[0];
		coord3[1] = coord_ptr[1];
		coord3[2] = coord_ptr[2];

		coord_ptr = a_uGrid->GetPoint(p4);
		coord4[0] = coord_ptr[0];
		coord4[1] = coord_ptr[1];
		coord4[2] = coord_ptr[2];

		bool bottomCoplane = fourPointsCoplanar(coord1,
			coord2,
			coord3,
			coord4);

		vtkIdType verts[2 * 4];
		//debug1<<"2D node for this hexahedron ";
		//first add bottom face node
		for (int iNode = 0; iNode < 4; ++iNode)
		{
			long p = a_validBottomNode[iNode];
			// debug1<<" "<<p;
			verts[iNode] = a_2DPointto3DPoints[p*a_num_layers + a_layerID];
		}
		//then add top face node
		for (int iNode = 4; iNode < 2 * 4; ++iNode)
		{
			//int p = a_validTopNode[iNode-4];
			// debug1<<" "<<p;
			//verts[iNode] = a_2DPointto3DPoints[p*m_num_layers+a_layerID+1];
			verts[iNode] = p_t[iNode - 4];
		}

		// this is a HEXAHEDRON /

		a_uGrid->InsertNextCell(VTK_HEXAHEDRON, 2 * 4, verts);
		//debug1<<"hexahedron at level (top and bottom coplane) " <<a_layerID<<" cell "<<a_Cell<<endl;
		m_quad_hexhedron++;
	}
	else if (a_validBottomNodeNum == 3) // still do it using hexahedron by using top point of degenerated point twice 
	{

		vtkIdType verts[8];

		//first add bottom face node
		for (int iNode = 0; iNode < 3; ++iNode)
		{
			long p = a_validBottomNode[iNode];
			// debug1<<" "<<p;
			verts[iNode] = a_2DPointto3DPoints[p*a_num_layers + a_layerID];
		}
		verts[3] = p_t[3];
		//then add top face node
		for (int iNode = 4; iNode < 8; ++iNode)
		{
			long p = a_validTopNode[iNode - 4];
			// debug1<<" "<<p;
			// verts[iNode] = a_2DPointto3DPoints[p*m_num_layers+a_layerID+1];
			verts[iNode] = p_t[iNode - 4];
		}
		a_uGrid->InsertNextCell(VTK_HEXAHEDRON, 2 * 4, verts);
		m_quad_hexhedron++;
	}
	else if (a_validBottomNodeNum == 2) // bottom have two node, this is a WEDGE
	{
		//debug1 <<"wedge at layer " <<a_layerID<<" cell "
		//		<<a_Cell<<" "<<topCoplane<<"\n";

		vtkIdType verts[6];
		int cellNodes[4];

		long p1 = a_validBottomNode[0];
		long p2 = a_validBottomNode[1];
		int p1_loc = 0;
		int p2_loc = 0;
		for (int i = 0; i < 4; i++)
		{
			cellNodes[i] = (a_faceNodePtr[i + 1] - 1);
			if (cellNodes[i] == p1)
			{
				p1_loc = i;
			}
			if (cellNodes[i] == p2)
			{
				p2_loc = i;
			}
		}



		// add first side of wedge
		verts[0] = a_2DPointto3DPoints[p1*a_num_layers + a_layerID];
		verts[1] = a_2DPointto3DPoints[p1*a_num_layers + a_layerID + 1];


		int p1NeiborFront = p1_loc + 1;
		if (p1NeiborFront >= 4)
		{
			p1NeiborFront = p1NeiborFront - 4;
		}
		int p1NeiborBack = p1_loc - 1;
		if (p1NeiborBack < 0)
		{
			p1NeiborBack = p1NeiborBack + 4;
		}
		int p3 = cellNodes[p1NeiborBack];
		if (p3 == p2)
		{
			p3 = cellNodes[p1NeiborFront];
		}

		int valid_bottom = std::max(1, a_kbp_node[p3]) - 1;
		if ((a_layerID + 1) < valid_bottom)
		{
			verts[2] = a_2DPointto3DPoints[p3*a_num_layers + valid_bottom];
		}
		else
		{
			verts[2] = a_2DPointto3DPoints[p3*a_num_layers + a_layerID + 1];
		}

		verts[3] = a_2DPointto3DPoints[p2*a_num_layers + a_layerID];
		verts[4] = a_2DPointto3DPoints[p2*a_num_layers + a_layerID + 1];

		int p2NeiborFront = p2_loc + 1;
		if (p2NeiborFront >= 4)
		{
			p2NeiborFront = p2NeiborFront - 4;
		}
		int p2NeiborBack = p2_loc - 1;
		if (p2NeiborBack < 0)
		{
			p2NeiborBack = p2NeiborBack + 4;
		}

		p3 = cellNodes[p2NeiborBack];
		if (p3 == p1)
		{
			p3 = cellNodes[p2NeiborFront];
		}

		valid_bottom = std::max(1, a_kbp_node[p3]) - 1;
		if ((a_layerID + 1) < valid_bottom)
		{
			verts[5] = a_2DPointto3DPoints[p3*a_num_layers + valid_bottom];
		}
		else
		{
			verts[5] = a_2DPointto3DPoints[p3*a_num_layers + a_layerID + 1];
		}


		a_uGrid->InsertNextCell(VTK_WEDGE, 6, verts);
		m_quad_wedge++;
	}
	else // this is pyramid 5 nodes
	{
		//debug1 <<"pyramid at layer " <<a_layerID<<" cell "
		//		<<a_Cell<<" "<<topCoplane<<"\n";;
		vtkIdType verts[5];

		for (int iNode = 0; iNode < 4; iNode++)
		{
			verts[iNode] = p_t[iNode];
		}
		long p3 = a_validBottomNode[0];
		verts[4] = a_2DPointto3DPoints[p3*a_num_layers + a_layerID];
		//debug1<<" bottom: "<<p3<<"\n";
		a_uGrid->InsertNextCell(VTK_PYRAMID, 5, verts);
		m_quad_pyramid++;
	}

}

void    avtMDSCHISMFileFormatImpl::validTopBottomNode(int       &   a_validTopNodeNum,
	int       &   a_validBottomNodeNum,
	long      *   a_validTopNode,
	long      *   a_validBottomNode,
	const int &   a_layerID,
	long      *   a_faceNodePtr,
	int       *   a_kbp_node ) const

{

	int numberOfNodeInCell = *a_faceNodePtr;


	//get node id indexed in 2d mesh
	long * node2D = new long[numberOfNodeInCell];

	for (int iNode = 0; iNode < numberOfNodeInCell; iNode++)
	{
		// faceNodePtr stores node index starting from 1
		// need decrease by one to be conisistent with
		// visit custom
		node2D[iNode] = a_faceNodePtr[iNode + 1] - 1;

	}


	//find out valid number node on bottom face
	int validBottomNode = 0;
	for (int iNode = 0; iNode < numberOfNodeInCell; iNode++)
	{
		// layer in m_kbp_node starts from 1, a_layerID starts from 0
		int bottomLayer = a_kbp_node[node2D[iNode]];
		if ((a_layerID + 1) >= bottomLayer)
		{
			a_validBottomNode[validBottomNode] = node2D[iNode];
			validBottomNode++;
		}
	}

	//find out valid number node on top face
	//int validTopNode = numberOfNodeInCell;
	//for(int iNode=0;iNode<numberOfNodeInCell; ++iNode)
	//{
	//   a_validTopNode[iNode]=node2D[iNode];

	//}
	//find out valid number node on top face
	long validTopNode = 0;
	for (int iNode = 0; iNode < numberOfNodeInCell; iNode++)
	{
		// layer in m_kbp_node starts from 1, a_layerID starts from 0
		//int bottomLayer = m_kbp_node[node2D[iNode]];
		//if ((a_layerID+2) >= bottomLayer)
		{
			a_validTopNode[validTopNode] = node2D[iNode];
			validTopNode++;
		}

	}
	a_validTopNodeNum = validTopNode;
	a_validBottomNodeNum = validBottomNode;
	delete node2D;

}
// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl:: updateMeshZCoordinates
//
//  Purpose:
//      update mesh points z coordinates for surface changes with time
//  
//  Arguments:
//     a_pointSet  the points making up mesh
//     a_timeState time
//     a_meshName  mesh label
//
//  Programmer: qshu
//  Creation:   Mon Sep 24 2012
//
// ****************************************************************************
void  avtMDSCHISMFileFormatImpl::updateMeshZCoordinates(vtkPoints * a_pointSet,
	const int   a_timeState,
	const int   a_domainid,
	const char* a_meshName)
{

	std::map<int, MDSCHISMMeshProvider*>::iterator it;
	it = m_external_mesh_providers.find(a_domainid);
	if (it == m_external_mesh_providers.end())
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve domain mesh provider of domain" << a_domainid;
		EXCEPTION1(InvalidVariableException, msgStream.str());
	}

	MDSCHISMMeshProvider * mesh_provider = it->second;

	if (!strcmp(a_meshName, m_mesh_2d.c_str()))
	{
		float * pointPtr = (float *)a_pointSet->GetVoidPointer(0);
		std::map<int, long>::iterator it;
		it = m_num_mesh_nodes.find(a_domainid);
		if (it == m_num_mesh_nodes.end())
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve num of 2d nodes of domain" << a_domainid;
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		long num_mesh_nodes = it->second;
		float * zCache = new float[num_mesh_nodes];
		if (!mesh_provider->zcoords2D(zCache, a_timeState))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve node z at step " << a_timeState << " domain " << a_domainid;
			EXCEPTION1(DBYieldedNoDataException, msgStream.str());
		}
		for (int iNode = 0; iNode < num_mesh_nodes; iNode++)
		{
			float z = zCache[iNode];
			pointPtr++;
			pointPtr++;
			*pointPtr++ = z;
			//debug1<<1<<" "<<iNode<<" "<<x<<" "<<y<<" "<<z<<"\n";
		}
		delete zCache;
	}
	else if (!strcmp(a_meshName, m_mesh_3d.c_str()))
	{
		float * pointPtr = (float *)a_pointSet->GetVoidPointer(0);
		std::map<int, long>::iterator it;
		it = m_total_valid_3D_point.find(a_domainid);
		if (it == m_total_valid_3D_point.end())
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve num of 3d nodes of domain" << a_domainid;
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		long total_valid_3D_point = it->second;

		float * nodeZPtr = new float[total_valid_3D_point];
		if (!mesh_provider->zcoords3D(nodeZPtr, a_timeState))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve node z at step " << a_timeState;
			EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
		}
		for (int iNode = 0; iNode < total_valid_3D_point; iNode++)
		{
			float z = nodeZPtr[iNode];
			pointPtr++;
			pointPtr++;
			*pointPtr++ = z;
		}
		delete nodeZPtr;
	}
	else if (!strcmp(a_meshName, m_side_center_point_2d_mesh.c_str()))
	{
		float * pointPtr = (float *)a_pointSet->GetVoidPointer(0);
		std::map<int, long>::iterator it;
		it = m_num_mesh_edges.find(a_domainid);
		if (it == m_num_mesh_edges.end())
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve num of 3d nodes of domain" << a_domainid;
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		long num_mesh_edges = it->second;
		float * zCache = new float[num_mesh_edges];
		if (!mesh_provider->zSideCenter2D(zCache, a_timeState))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve node z at step " << a_timeState;
			EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
		}
		for (int i = 0; i < num_mesh_edges; i++)
		{
			float z = zCache[i];
			pointPtr++;
			pointPtr++;
			*pointPtr++ = z;
			//debug1<<1<<" "<<iNode<<" "<<x<<" "<<y<<" "<<z<<"\n";
		}
		delete zCache;
	}
	else if (!strcmp(a_meshName, m_side_center_point_3d_mesh.c_str()))
	{
		float * pointPtr = (float *)a_pointSet->GetVoidPointer(0);
		std::map<int, long>::iterator it;
		it = m_total_valid_3D_side.find(a_domainid);
		if (it == m_total_valid_3D_side.end())
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve num of 3d side of domain" << a_domainid;
			EXCEPTION1(InvalidVariableException, msgStream.str());
		}
		long total_valid_3D_side = it->second;
		float * sideCenterZPtr = new float[total_valid_3D_side];
		//loadAndCacheZSide(a_timeState,sideCenterZPtr);
		if (!mesh_provider->zSideCenter3D(sideCenterZPtr, a_timeState))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve side z at step " << a_timeState;
			EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
		}
		for (int i = 0; i < total_valid_3D_side; i++)
		{
			float z = sideCenterZPtr[i];
			pointPtr++;
			pointPtr++;
			*pointPtr++ = z;
		}
		delete sideCenterZPtr;
	}



}

void avtMDSCHISMFileFormatImpl::loadAndCacheZSide(const int& a_timeState, float * a_sideCenterZPtr)
{
	return;
}

void avtMDSCHISMFileFormatImpl::loadAndCacheZEle(const int& a_timeState, float * a_eleCenterZPtr)
{
	return;

}

void avtMDSCHISMFileFormatImpl::loadAndCacheZ(const int& a_timeState, float* a_nodeZPtr)
{
	return;
}


void avtMDSCHISMFileFormatImpl::load_node_dry_wet(const int & a_time, const int & a_domain)
{

	long num = m_num_mesh_nodes[a_domain];

	std::map<int, int>::iterator it;
	it = m_node_dry_wet_cached_time.find(a_domain);
	int node_dry_wet_cached_time = MeshConstants10::INVALID_NUM;

	if (it != (m_node_dry_wet_cached_time.end()))
	{
		node_dry_wet_cached_time = it->second;
	}
	std::map<int, int*>::iterator it2;
		it2 = m_node_dry_wet.find(a_domain);
		int * node_dry_wet = NULL;
		if (it2 != (m_node_dry_wet.end()))
		{
			node_dry_wet = it2->second;
		}

	if ((node_dry_wet_cached_time == a_time) && (node_dry_wet))
	{
		return;
	}
	else
	{
		if (!(node_dry_wet))
		{
			node_dry_wet = new int[num];
			m_node_dry_wet[a_domain] = node_dry_wet;
		}

		//make sure ele dry wet is loaded and cached for this time before fill nodes dry/wet
		std::map<int, int*>::iterator it3;
		it3 = m_ele_dry_wet.find(a_domain);
		int * ele_dry_wet = NULL;
		if (it3 != (m_ele_dry_wet.end()))
		{
			ele_dry_wet = it3->second;
		}

		if (!ele_dry_wet)
		{
			EXCEPTION1(DBYieldedNoDataException, "fill node dry/wet before fill element dry/wet \n");
		}

		m_external_mesh_providers[a_domain]->fill_node_dry_wet(node_dry_wet, ele_dry_wet);

			long nominal_data_size_per_layer = 0;


		for (int i = 0; i < num; i++)
		{
			if (!(node_dry_wet[i]))
			{
				nominal_data_size_per_layer++;
			}
		}

		m_nominal_size_per_layer[MeshConstants10::NODE + std::to_string(a_domain)] = nominal_data_size_per_layer;
		m_node_dry_wet_cached_time[a_domain] = a_time;
	}

	return;
}


void avtMDSCHISMFileFormatImpl::load_ele_dry_wet(const int & a_time, const int& a_domain)
{
	long num = m_num_mesh_faces[a_domain];
	std::map<int, int>::iterator it;
	it = m_ele_dry_wet_cached_time.find(a_domain);
	int ele_dry_wet_cached_time = MeshConstants10::INVALID_NUM;

	if (it != (m_ele_dry_wet_cached_time.end()))
	{
		ele_dry_wet_cached_time = it->second;
	}

	if (ele_dry_wet_cached_time == a_time)
	{
		return;
	}
	else
	{
		int * ele_dry_wet = NULL;
		std::map<int, int*>::iterator it3;
		it3 = m_ele_dry_wet.find(a_domain);

		if (it3 != (m_ele_dry_wet.end()))
		{
			ele_dry_wet = it3->second;
		}
		if (!ele_dry_wet)
		{
			ele_dry_wet = new int[num];
			m_ele_dry_wet[a_domain] = ele_dry_wet;
		}

		//debug1 << "done allocate for ele dry/wet for domain " << a_domain <<" "<<num<< "\n";
		std::string SCHISMVarName = MeshConstants10::ELEM_DRYWET;
		SCHISMVar10* SCHISMVarPtr = NULL;
		MDSchismOutput * data_file_ptr = NULL;

		std::map<int, MDSchismOutput*>::iterator it4;
		it4 = m_data_files.find(a_domain);

		if (it4 != (m_data_files.end()))
		{
			data_file_ptr = it4->second;
		}
		if (!data_file_ptr)
		{
			EXCEPTION1(InvalidVariableException, ("invalid data file " + m_data_file).c_str());
		}

		//debug1 << "begin load ele dry/wet for domain " << a_domain << "\n";
		try
		{
			SCHISMVarPtr = data_file_ptr->get_var(SCHISMVarName);
		}
		catch (...)
		{
			EXCEPTION1(InvalidVariableException, ("no " + SCHISMVarName + " in " + m_data_file).c_str());
		}
		//debug1 << "get var ele dry/wet for domain " << a_domain << "\n";
		SCHISMVarPtr->set_cur(a_time);
		
		//original dry wet is defined as float, use a float buffer here
		float * buffer = new float[num];

		if (!(SCHISMVarPtr->get(buffer)))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve " << SCHISMVarName << " at step " << a_time << " domain " << a_domain;
			EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
		}
		for (int i = 0; i < num; i++)
		{
			ele_dry_wet[i] = int(buffer[i]);
			
		}
		delete buffer;
		long nominal_data_size_per_layer = 0;
		//debug1 << "done load ele dry/wet for domain " << a_domain << "\n";
		for (int i = 0; i < num; i++)
		{
			if (!(ele_dry_wet[i]))
			{
				nominal_data_size_per_layer++;
			}
		}

		m_nominal_size_per_layer[MeshConstants10::ELEM + std::to_string(a_domain)] = nominal_data_size_per_layer;
		m_ele_dry_wet_cached_time[a_domain] = a_time;
	}

	return;
}


void avtMDSCHISMFileFormatImpl::load_side_dry_wet(const int & a_time, const int& a_domain)
{
	long num = m_num_mesh_edges[a_domain];

	std::map<int, int>::iterator it;
	it = m_side_dry_wet_cached_time.find(a_domain);
	int side_dry_wet_cached_time = MeshConstants10::INVALID_NUM;

		if (it != (m_side_dry_wet_cached_time.end()))
		{
			side_dry_wet_cached_time = it->second;
		}

	std::map<int, int*>::iterator it2;
		it2 = m_side_dry_wet.find(a_domain);
		int * side_dry_wet = NULL;
		if (it2 != (m_side_dry_wet.end()))
		{
			side_dry_wet = it2->second;
		}

	if ((side_dry_wet_cached_time == a_time) && (side_dry_wet))
	{
		return;
	}
	else
	{
		int * ele_dry_wet = NULL;
		std::map<int, int*>::iterator it3;
		it3 = m_ele_dry_wet.find(a_domain);

		if (it3 != (m_ele_dry_wet.end()))
		{
			ele_dry_wet = it3->second;
		}

		if (!ele_dry_wet)
		{
			EXCEPTION1(DBYieldedNoDataException, "fill side dry/wet before fill element dry/wet \n");
		}

		if (!(side_dry_wet))
		{
			side_dry_wet = new int[num];
			m_side_dry_wet[a_domain] = side_dry_wet;
		}


		m_external_mesh_providers[a_domain]->fill_side_dry_wet(side_dry_wet, ele_dry_wet);

		long nominal_data_size_per_layer = 0;

		for (int i = 0; i < num; i++)
		{
			if (!(side_dry_wet[i]))
			{
				nominal_data_size_per_layer++;
			}
		}

		m_nominal_size_per_layer[MeshConstants10::EDGE + std::to_string(a_domain)] = nominal_data_size_per_layer;

		m_side_dry_wet_cached_time[a_domain] = a_time;
	}

	return;
}

void   avtMDSCHISMFileFormatImpl::getMeshDimensions(const int& a_domainID, MeshProvider10 * a_meshProviderPtr)
{

	//debug1 << "getting mesh node number\n";
	long num_mesh_nodes = MeshConstants10::INVALID_NUM;
	num_mesh_nodes = a_meshProviderPtr->numberOfNode();

	//debug1 << "number of node " << m_num_mesh_nodes[a_domainID] << "\n";

	if (num_mesh_nodes < 0)  //to do: face centered data will get num faces
	{
		EXCEPTION1(InvalidVariableException, ("no nodes found in " + a_meshProviderPtr->file()).c_str());
	}
	m_num_mesh_nodes[a_domainID] = num_mesh_nodes;

	// SCHISMDim constructor and destructor are private, no need to del.
	long num_mesh_faces = MeshConstants10::INVALID_NUM;
	num_mesh_faces = a_meshProviderPtr->numberOfElement();

	if (num_mesh_faces < 0)
	{
		EXCEPTION1(InvalidVariableException, ("no element found in " + a_meshProviderPtr->file()).c_str());
	}
	m_num_mesh_faces[a_domainID] = num_mesh_faces;
	//debug1 << "get face num" << num_mesh_faces << "\n";

	long num_mesh_edges = MeshConstants10::INVALID_NUM;
	num_mesh_edges = a_meshProviderPtr->numberOfSide();

	if (num_mesh_edges < 0)
	{
		EXCEPTION1(InvalidVariableException, ("no edges found in " + a_meshProviderPtr->file()).c_str());
	}
	m_num_mesh_edges[a_domainID] = num_mesh_edges;
	//initlaize nominal data size per layer for all centering
	m_nominal_size_per_layer[MeshConstants10::NODE + std::to_string(a_domainID)] = m_num_mesh_nodes[a_domainID];
	m_nominal_size_per_layer[MeshConstants10::ELEM + std::to_string(a_domainID)] = m_num_mesh_faces[a_domainID];
	m_nominal_size_per_layer[MeshConstants10::EDGE + std::to_string(a_domainID)] = m_num_mesh_edges[a_domainID];

	m_num_layers = a_meshProviderPtr->numberOfLayer();

	if (m_num_layers < 0)
	{
		EXCEPTION1(InvalidVariableException, ("no layer found ,one required at least, in " + a_meshProviderPtr->file()).c_str());
	}
	//debug1 << "get layers num" << m_num_layers << "\n";
	debug1 << "done load mesh dimension for mesh block " << a_domainID << "\n";
}

void    avtMDSCHISMFileFormatImpl::load_bottom(const int& a_time, const int & a_domain)
{
	int cache_node_kbp_id = MeshConstants10::INVALID_NUM;

	long num = m_num_mesh_nodes[a_domain];

	std::map<int, int>::iterator it;
	it = m_kbp_node_time.find(a_domain);
	std::map<int, int*>::iterator it2;
	it2 = m_kbp_node.find(a_domain);

	if (it != m_kbp_node_time.end())
	{
		cache_node_kbp_id = it->second;
	}
	int * kbp_node = NULL;
	if (it2 != m_kbp_node.end())
	{
		kbp_node = it2->second;
	}


	if ((m_cache_kbp_id != a_time) || (!kbp_node))
	{
		if (!kbp_node)
		{
			kbp_node = new int[num];
			m_kbp_node[a_domain] = kbp_node;
		}

		m_external_mesh_providers[a_domain]->fillKbp00(kbp_node, a_time);
		m_kbp_node_time[a_domain] = a_time;
	}


	//count total valid 3d point
    std:map < long, int> domain_3d_to_2d;
	long total_valid_3D_point = 0;
	for (int iLayer = 0; iLayer < m_num_layers; iLayer++)
	{
		for (int iNode = 0; iNode < m_num_mesh_nodes[a_domain]; iNode++)
		{
			int bottomLayer = kbp_node[iNode];
			if (bottomLayer <= (iLayer + 1))
			{
				domain_3d_to_2d[total_valid_3D_point] = iNode;
				total_valid_3D_point++;
			}
		}
	}

	m_total_valid_3D_point[a_domain] = total_valid_3D_point;
	m_3d_node_to_2d_node[a_domain] = domain_3d_to_2d;


	int cache_ele_kbp_id = MeshConstants10::INVALID_NUM;

	num = m_num_mesh_faces[a_domain];
	it = m_kbp_ele_time.find(a_domain);
	it2 = m_kbp_ele.find(a_domain);

	if (it != m_kbp_ele_time.end())
	{
		cache_ele_kbp_id = it->second;
	}
	int * kbp_ele = NULL;
	if (it2 != m_kbp_ele.end())
	{
		kbp_ele = it2->second;
	}


	if ((cache_ele_kbp_id != a_time) || (!kbp_ele))
	{
		if (!kbp_ele)
		{
			kbp_ele = new int[num];
			m_kbp_ele[a_domain] = kbp_ele;
		}

		m_external_mesh_providers[a_domain]->fillKbe(kbp_ele, a_time);
		m_kbp_ele_time[a_domain] = a_time;
	}


	//count total valid 3d ele
	long total_valid_3D_ele = 0;
	for (int iLayer = 0; iLayer < m_num_layers; iLayer++)
	{
		for (int iNode = 0; iNode < m_num_mesh_faces[a_domain]; iNode++)
		{
			int bottomLayer = kbp_ele[iNode];
			if (bottomLayer <= (iLayer + 1))
			{
				total_valid_3D_ele++;
			}
		}
	}

	m_total_valid_3D_ele[a_domain] = total_valid_3D_ele;

	int cache_side_kbp_id = MeshConstants10::INVALID_NUM;
	num = m_num_mesh_edges[a_domain];
	it = m_kbp_side_time.find(a_domain);
	it2 = m_kbp_side.find(a_domain);

	if (it != m_kbp_side_time.end())
	{
		cache_side_kbp_id = it->second;
	}
	int * kbp_side = NULL;
	if (it2 != m_kbp_side.end())
	{
		kbp_side = it2->second;
	}


	if ((cache_side_kbp_id != a_time) || (!kbp_side))
	{
		if (!kbp_side)
		{
			kbp_side = new int[num];
			m_kbp_side[a_domain] = kbp_side;
		}

		m_external_mesh_providers[a_domain]->fillKbs(kbp_side, a_time);
		m_kbp_side_time[a_domain] = a_time;
	}


	//count total valid 3d ele
	long total_valid_3D_side = 0;
	for (int iLayer = 0; iLayer < m_num_layers; iLayer++)
	{
		for (int iNode = 0; iNode < m_num_mesh_edges[a_domain]; iNode++)
		{
			int bottomLayer = kbp_side[iNode];
			if (bottomLayer <= (iLayer + 1))
			{
				total_valid_3D_side++;
			}
		}
	}

	m_total_valid_3D_side[a_domain] = total_valid_3D_side;

}



// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::getLayer
//
//  Purpose:
//      Return layer id for each 3d prsim
//      
//     
//
//  Arguments:
//       
//                  
//      None
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Mon DEC 13 09:57:49 PDT 2015
//
// ****************************************************************************

vtkDataArray*   avtMDSCHISMFileFormatImpl::getLayer(const int& a_domain)
{
	int * kbp00 = m_kbp_ele[a_domain];
	std::string data_center = MeshConstants10::ELEM;

	int valid_var_size = 0;
	int num_total_layers = m_num_layers - 1;

	long nominal_num_data_per_Layer = m_num_mesh_faces[a_domain];

	for (long iNode = 0; iNode < nominal_num_data_per_Layer; iNode++)
	{
		//if(!(m_ele_dry_wet[iNode]))
		{
			valid_var_size += num_total_layers - std::max(1, kbp00[iNode]) + 1;
		}
	}

	long ntuples = valid_var_size;
	vtkIntArray *rv = vtkIntArray::New();
	rv->SetNumberOfTuples(ntuples);
	long idata = 0;


	for (int iLayer = 0; iLayer < num_total_layers; iLayer++)
	{
		for (long iEle = 0; iEle < nominal_num_data_per_Layer; iEle++)
		{
			int valid_bottom_layer = std::max(1, kbp00[iEle]) - 1;
			if (iLayer >= valid_bottom_layer)
			{
				rv->SetTuple1(idata, iLayer);
				idata++;
			}
		}
	}

	return rv;
}

// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::getLayer
//
//  Purpose:
//      Return layer id for 3d node, ele or side according to name
//      
//     
//
//  Arguments:
//       
//                  
//      level id in string
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Mon DEC 13 09:57:49 PDT 2015
//
// ****************************************************************************

vtkDataArray*   avtMDSCHISMFileFormatImpl::getLevel(const string& a_level_name, const int& a_domain)
{
	int * kbp00;
	std::string data_center(MeshConstants10::NODE);
	long nominal_num_data_per_Layer = 0;
	if (!(a_level_name.compare(MeshConstants10::NODE_LEVEL)))
	{
		//kbp00=m_kbp_data;
		kbp00 = m_kbp_node[a_domain];
		nominal_num_data_per_Layer = m_num_mesh_nodes[a_domain];
	}
	else if (!(a_level_name.compare(MeshConstants10::ELE_LEVEL)))
	{
		kbp00 = m_kbp_ele[a_domain];
		data_center = MeshConstants10::ELEM;
		nominal_num_data_per_Layer = m_num_mesh_faces[a_domain];
	}
	else if (!(a_level_name.compare(MeshConstants10::SIDE_LEVEL)))
	{
		kbp00 = m_kbp_side[a_domain];
		data_center = MeshConstants10::EDGE;
		nominal_num_data_per_Layer = m_num_mesh_edges[a_domain];
	}
	else
	{
		stringstream msgStream(stringstream::out);
		msgStream << a_level_name << " is tot a valid mesh level variable\n";
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}

	int valid_var_size = 0;
	int num_total_layers = m_num_layers;

	for (long iNode = 0; iNode < nominal_num_data_per_Layer; iNode++)
	{

		valid_var_size += num_total_layers - std::max(1, kbp00[iNode]) + 1;
	}

	long ntuples = valid_var_size;
	vtkIntArray *rv = vtkIntArray::New();
	rv->SetNumberOfTuples(ntuples);
	long idata = 0;


	for (int iLayer = 0; iLayer < num_total_layers; iLayer++)
	{
		for (long iNode = 0; iNode < nominal_num_data_per_Layer; iNode++)
		{
			int valid_bottom_layer = std::max(1, kbp00[iNode]) - 1;
			if (iLayer >= valid_bottom_layer)
			{
				rv->SetTuple1(idata, iLayer);
				idata++;
			}
		}
	}

	return rv;
}

vtkDataArray*   avtMDSCHISMFileFormatImpl::getBottom(const string& a_bottom_name, const int& a_domain)
{
	int * kbp00;
	std::string data_center(MeshConstants10::NODE);
	long nominal_num_data_per_Layer = 0;
	if (!(a_bottom_name.compare(MeshConstants10::NODE_BOTTOM)))
	{
		//kbp00=m_kbp_data;
		kbp00 = m_kbp_node[a_domain];
		nominal_num_data_per_Layer = m_num_mesh_nodes[a_domain];
	}
	else if (!(a_bottom_name.compare(MeshConstants10::FACE_BOTTOM)))
	{
		kbp00 = m_kbp_ele[a_domain];
		data_center = MeshConstants10::ELEM;
		nominal_num_data_per_Layer = m_num_mesh_faces[a_domain];
	}
	else if (!(a_bottom_name.compare(MeshConstants10::EDGE_BOTTOM)))
	{
		kbp00 = m_kbp_side[a_domain];
		data_center = MeshConstants10::EDGE;
		nominal_num_data_per_Layer = m_num_mesh_edges[a_domain];
	}
	else
	{
		stringstream msgStream(stringstream::out);
		msgStream << a_bottom_name << " is not a valid mesh bottom variable\n";
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}


	//long nominal_num_data_per_Layer=m_nominal_size_per_layer[data_center];


	long ntuples = nominal_num_data_per_Layer;
	vtkIntArray *rv = vtkIntArray::New();
	rv->SetNumberOfTuples(ntuples);
	long idata = 0;



	for (long iNode = 0; iNode < nominal_num_data_per_Layer; iNode++)
	{

		rv->SetTuple1(idata, kbp00[iNode]);
		idata++;

	}


	return rv;
}

vtkDataArray*   avtMDSCHISMFileFormatImpl::get_node_global_id( const int& a_domain)
{
	
	std::string data_center(MeshConstants10::NODE);
	long nominal_num_data_per_layer = 0;
	nominal_num_data_per_layer = m_num_mesh_nodes[a_domain];


	long ntuples = nominal_num_data_per_layer;
	vtkIntArray *rv = vtkIntArray::New();
	rv->SetNumberOfTuples(ntuples);
	long idata = 0;
	long * buff = new long[nominal_num_data_per_layer];
	m_external_mesh_providers[a_domain]->fill_node_global_id(buff);

	for (long iNode = 0; iNode < nominal_num_data_per_layer; iNode++)
	{

		rv->SetTuple1(idata, buff[iNode]);
		idata++;

	}


	return rv;
}


vtkDataArray*   avtMDSCHISMFileFormatImpl::get_ele_global_id(const int& a_domain)
{
	
	std::string data_center(MeshConstants10::ELEM);

	long nominal_num_data_per_layer = 0;
	
    data_center = MeshConstants10::ELEM;
    nominal_num_data_per_layer = m_num_mesh_faces[a_domain];

	long ntuples = nominal_num_data_per_layer;
	vtkIntArray *rv = vtkIntArray::New();
	rv->SetNumberOfTuples(ntuples);
	long idata = 0;
	long * buff = new long[nominal_num_data_per_layer];
	m_external_mesh_providers[a_domain]->fill_ele_global_id(buff);
	for (long iEle = 0; iEle < nominal_num_data_per_layer; iEle++)
	{
		rv->SetTuple1(idata,buff[iEle]);
		idata++;

	}
	delete buff;

	return rv;
}
vtkDataArray*   avtMDSCHISMFileFormatImpl::get_node_depth(const int& a_domain)
{

	std::string data_center(MeshConstants10::NODE);
	long nominal_num_data_per_layer = 0;
	nominal_num_data_per_layer = m_num_mesh_nodes[a_domain];


	long ntuples = nominal_num_data_per_layer;
	vtkDoubleArray *rv = vtkDoubleArray::New();
	rv->SetNumberOfTuples(ntuples);
	long idata = 0;
	double * buff = new double[nominal_num_data_per_layer];
	m_external_mesh_providers[a_domain]->depth(buff);

	for (long iNode = 0; iNode < nominal_num_data_per_layer; iNode++)
	{

		rv->SetTuple1(idata, buff[iNode]);
		idata++;

	}


	return rv;
}
// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      varname    The name of the variable requested.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Wed Mar 13 09:13:49 PDT 2013
//
// ****************************************************************************

vtkDataArray *
avtMDSCHISMFileFormatImpl::GetVar(int a_timeState, int a_domainID, const char *a_varName)
{

	//a_varName is the label used by visit. This name
 //needs to be map into varname used in data file.

	if (!strcmp(a_varName, (MeshConstants10::LAYER).c_str()))
	{
		return getLayer(a_domainID);
	}
	else  if (!strcmp(a_varName, (MeshConstants10::NODE_LEVEL).c_str()))
	{
		load_node_dry_wet(a_timeState, a_domainID);
		return getLevel(MeshConstants10::NODE_LEVEL, a_domainID);
	}
	else  if (!strcmp(a_varName, (MeshConstants10::ELE_LEVEL).c_str()))
	{
		load_ele_dry_wet(a_timeState, a_domainID);
		return getLevel(MeshConstants10::ELE_LEVEL, a_domainID);
	}
	else  if (!strcmp(a_varName, (MeshConstants10::SIDE_LEVEL).c_str()))
	{
		load_side_dry_wet(a_timeState, a_domainID);
		return getLevel(MeshConstants10::SIDE_LEVEL, a_domainID);
	}
	else  if (!strcmp(a_varName, (MeshConstants10::NODE_BOTTOM).c_str()))
	{
		return getBottom(MeshConstants10::NODE_BOTTOM, a_domainID);
	}
	else  if (!strcmp(a_varName, (MeshConstants10::FACE_BOTTOM).c_str()))
	{
		return getBottom(MeshConstants10::FACE_BOTTOM, a_domainID);
	}
	else  if (!strcmp(a_varName, (MeshConstants10::EDGE_BOTTOM).c_str()))
	{
		return getBottom(MeshConstants10::EDGE_BOTTOM, a_domainID);
	}
	else  if (!strcmp(a_varName, (MeshConstants10::NODE_GLOBAL_ID).c_str()))
	{
		return get_node_global_id(a_domainID);
	}
	else  if (!strcmp(a_varName, (MeshConstants10::ELE_GLOBAL_ID).c_str()))
	{
		return get_ele_global_id(a_domainID);
	}
	else  if (!strcmp(a_varName, (MeshConstants10::NODE_DEPTH).c_str()))
	{
		return get_node_depth(a_domainID);
	}


	std::string SCHISMVarName = m_var_name_label_map[a_varName];
	std::string varMesh = m_var_mesh_map[a_varName];

	SCHISMVar10 * SCHISMVarPtr = m_data_files[a_domainID]->get_var(SCHISMVarName);

	std::string level_center = SCHISMVarPtr->get_vertical_center();
	std::string data_center = SCHISMVarPtr->get_horizontal_center();

	int * drywet;
	debug1 << "getting " << a_varName << " " << level_center << " " << data_center << "\n";

	long numDataPerLayer = m_num_mesh_nodes[a_domainID];
	if (data_center == MeshConstants10::ELEM)
	{
		load_ele_dry_wet(a_timeState, a_domainID);
		numDataPerLayer = m_num_mesh_faces[a_domainID];
		drywet = m_ele_dry_wet[a_domainID];
	}
	else if (data_center == MeshConstants10::EDGE)
	{
		load_side_dry_wet(a_timeState, a_domainID);
		numDataPerLayer = m_num_mesh_edges[a_domainID];
		drywet = m_side_dry_wet[a_domainID];
	}
	else
	{
		debug1 << "try to load node dry wet\n";
		load_node_dry_wet(a_timeState, a_domainID);
		drywet = m_node_dry_wet[a_domainID];
	}


	long record_size_adjust = 0;
	int layer_size_adjust = 0;
	bool is_half_layer = false;

	if (level_center == MeshConstants10::HALF_LAYER)
	{
		record_size_adjust = -1;
		layer_size_adjust = -1;
		is_half_layer = true;
	}

	if (!(SCHISMVarPtr->is_valid()))
	{

		EXCEPTION1(InvalidVariableException, a_varName);
	}




	if ((!strcmp(a_varName, m_node_depth_label.c_str())) ||
		(!SCHISMVarIs3D(SCHISMVarPtr)))
	{
		float * valBuff;
		int numData = numDataPerLayer;

		valBuff = new float[numData];

		getSingleLayerVar(valBuff,
			m_data_files[a_domainID],
			a_timeState,
			SCHISMVarName);


		//total number of data = nodes for a time step
		int ntuples = numData;
		vtkDoubleArray *rv = vtkDoubleArray::New();
		rv->SetNumberOfTuples(ntuples);
		int idata = 0;
		for (int iNode = 0; iNode < numData; iNode++)
		{
			float valTemp = valBuff[iNode];
			if ((!(drywet[iNode])) || (!strcmp(a_varName, m_node_depth_label.c_str())))
			{

				rv->SetTuple1(idata, valTemp);
			}
			else
			{
				if (m_dry_wet_flag)
					rv->SetTuple1(idata, MeshConstants10::DRY_STATE);
				else
					rv->SetTuple1(idata, valTemp);
			}

			idata++;
		}
		delete   valBuff;
		return rv;
	}



	debug1 << "get SCHISM var " << SCHISMVarName << "\n";

	int num_data_layers = m_num_layers + layer_size_adjust;

	std::string        varName(a_varName);
	int * layerStarts;
	bool is_bottom = false;
	bool is_surface = false;
	if (varName.find(m_surface_state_suffix) != string::npos)
	{
		num_data_layers = 1;
		layerStarts = new int[num_data_layers];
		layerStarts[0] = m_num_layers - 1 + layer_size_adjust;
		is_surface = true;
	}
	else if (varName.find(m_bottom_state_suffix) != string::npos)
	{
		num_data_layers = 1;
		layerStarts = new int[num_data_layers];
		//this is the layer next to the bottom assumed,it will change
		//for each node later
		layerStarts[0] = 1;
		is_bottom = true;
	}
	else
	{
		layerStarts = new int[num_data_layers];
		for (int iLayers = 0; iLayers < num_data_layers; iLayers++)
		{
			layerStarts[iLayers] = iLayers;
		}
	}


	float * valBuff;
	int numTotalLayers = m_num_layers+ layer_size_adjust;
	valBuff = new float[numTotalLayers*numDataPerLayer];


	int numOfRecord = 1;
	int nodeStart = 0;
	int timeStart = a_timeState;

	SCHISMVarPtr->set_cur(timeStart);

	debug1 << "loading data for var " << a_varName << " for domain " << a_domainID <<" "<< numDataPerLayer << "\n";

	if (!(SCHISMVarPtr->get(valBuff)))
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve " << a_varName << " at step " << a_timeState;
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}
	
	long valid_var_size = 0;
	long * node_start_index = new long[numDataPerLayer];
	long * num_data_at_layer = new long[numTotalLayers];

	for (int iLayer = 0; iLayer < numTotalLayers; iLayer++)
	{
		num_data_at_layer[iLayer] = 0;
	}

	int * kbp00 = m_kbp_node[a_domainID];

	if (!(data_center.compare(MeshConstants10::EDGE)))
	{
		kbp00 = m_kbp_side[a_domainID];
	}
	else if (!(data_center.compare(MeshConstants10::ELEM)))
	{
		if (is_half_layer)
		{
			kbp00 = m_kbp_prism[a_domainID];
		}
		else
		{
			kbp00 = m_kbp_ele[a_domainID];
		}
	}
	debug1 << "kbp loaded for domain " << a_domainID << "\n";

	int num_layers = m_num_layers + layer_size_adjust;

	long half_remove_cell = 0;
	if (is_half_layer) half_remove_cell = 1;
	for (int iNode = 0; iNode < numDataPerLayer; iNode++)
	{

		node_start_index[iNode] = valid_var_size;
		//if (!(drywet[iNode]))
		{
			long record_len = num_layers - std::max(1, kbp00[iNode]) + 1;
			
			//if((drywet[iNode])&&(record_len>0)) half_remove_cell+=record_len;
			//if(record_len<0) record_len=0; // only happened for dry element
			valid_var_size += record_len;

		}
	}


	vtkFloatArray *rv = vtkFloatArray::New();
	int idata = 0;
	if (varName.find(m_depth_average_suffix) == string::npos)
	{
		idata = 0;

		// count total number of data 
		int ntuples = 0;

		ntuples = valid_var_size;
		if (is_bottom || is_surface) //for surface and bottom data, same as number of 2D node
		{
			ntuples = numDataPerLayer;
		}
		//if (is_half_layer) ntuples -= half_remove_cell;
		rv->SetNumberOfTuples(ntuples);
		debug1 << " count total num of data on domain "<<a_domainID<<"  is " << ntuples << "\n ";
		for (int iLayer = 0; iLayer < num_data_layers; iLayer++)
		{
			int layer = layerStarts[iLayer];
			for (int iNode = 0; iNode < numDataPerLayer; iNode++)
			{
				int valid_bottom_layer = std::max(1, kbp00[iNode]) - 1;
				//if (is_half_layer) valid_bottom_layer++;
				if (is_bottom)
				{
					layer = valid_bottom_layer + 1;
					if (layer > (numTotalLayers - 1))
					{
						layer = numTotalLayers - 1;
					}
				}
				float valTemp = MeshConstants10::DUMMY_ELEVATION;
				int start_index = node_start_index[iNode];

				if (layer >= valid_bottom_layer)
				{
					long data_loc = start_index + layer - valid_bottom_layer;
					valTemp = valBuff[data_loc];
					if (!(drywet[iNode]))
					{

						rv->SetTuple1(idata, valTemp);
					}
					else
					{
						if (m_dry_wet_flag)
						{
							rv->SetTuple1(idata, MeshConstants10::DRY_STATE);
						}
						else
						{
							rv->SetTuple1(idata, valTemp);
						}
					}
					idata++;
				}

			}
		}

	}
	else
	{
		idata = 0;
		//total number of data =nodes for a time step
		int ntuples = numDataPerLayer;
		rv->SetNumberOfTuples(ntuples);

		float * averageState = new float[numDataPerLayer];
		debug1 << "begin averaging state\n";
		depthAverage(averageState, valBuff, node_start_index, a_timeState, a_domainID, data_center, level_center);
		debug1 << "total num of averaged data " << numDataPerLayer << "\n";
		idata = 0;
		for (int iNode = 0; iNode < numDataPerLayer; iNode++)
		{
			//if (!(drywet[iNode]))
			if ((!drywet[iNode]) || (drywet[iNode] && (!(m_dry_wet_flag))))
			{
				float valTemp = averageState[iNode];
				rv->SetTuple1(idata, valTemp);
			}
			else
			{
				rv->SetTuple1(idata, MeshConstants10::DRY_STATE);
			}

			idata++;
		}
		delete averageState;
	}

	delete valBuff;
	delete layerStarts;
	delete num_data_at_layer;
	debug1 << "done load data for " << SCHISMVarName << "\n";
	return rv;

}


// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      varname    The name of the variable requested.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Wed Mar 13 09:13:49 PDT 2013
//
// ****************************************************************************

vtkDataArray *
avtMDSCHISMFileFormatImpl::GetVectorVar(int a_timeState, int a_domainID, const char *a_varName)
{

	std::string SCHISMVarName = m_var_name_label_map[a_varName];
	std::string varMesh = m_var_mesh_map[a_varName];

	SCHISMVar10 * SCHISMVarPtr = m_data_files[a_domainID]->get_var(SCHISMVarName);
	if (!(SCHISMVarPtr->is_valid()))
	{

		EXCEPTION1(InvalidVariableException, a_varName);
	}
	std::string level_center = SCHISMVarPtr->get_vertical_center();
	std::string data_center = SCHISMVarPtr->get_horizontal_center();

	int * drywet;

	long numDataPerLayer = m_num_mesh_nodes[a_domainID];
	if (data_center == MeshConstants10::ELEM)
	{
		load_ele_dry_wet(a_timeState, a_domainID);
		numDataPerLayer = m_num_mesh_faces[a_domainID];
		drywet = m_ele_dry_wet[a_domainID];
	}
	else if (data_center == MeshConstants10::EDGE)
	{
		load_side_dry_wet(a_timeState, a_domainID);
		numDataPerLayer = m_num_mesh_edges[a_domainID];
		drywet = m_side_dry_wet[a_domainID];
	}
	else
	{
		debug1 << "try to load node dry wet\n";
		load_node_dry_wet(a_timeState, a_domainID);
		drywet = m_node_dry_wet[a_domainID];
	}

	// last dim is vector component  
	int      numDim = SCHISMVarPtr->num_dims();
	SCHISMDim10* comDim = SCHISMVarPtr->get_dim(numDim - 1);
	int ncomps = comDim->size();
	int ucomps = (ncomps == 2 ? 3 : ncomps);
	float *oneEntry = new float[ucomps];
	int idata = 0;

	if (!SCHISMVarIs3D(SCHISMVarPtr))
	{
		float * valBuff;
		long numData = numDataPerLayer;

		long numDataEntry = numData;

		numData *= ncomps;

		valBuff = new float[numData];
		getSingleLayerVar(valBuff,
			m_data_files[a_domainID],
			a_timeState,
			SCHISMVarName);
		//total number of data = nodes for a time step
		long ntuples = numDataEntry;
		vtkFloatArray *rv = vtkFloatArray::New();
		rv->SetNumberOfComponents(ucomps);
		rv->SetNumberOfTuples(ntuples);
		idata = 0;
		for (long iNode = 0; iNode < ntuples; iNode++)
		{

			if ((!drywet[iNode]) || (drywet[iNode] && (!(m_dry_wet_flag))))
			{
				for (int iComp = 0; iComp < ncomps; iComp++)
				{
					oneEntry[iComp] = valBuff[iNode*ncomps + iComp];
				}
				for (int iComp = ncomps; iComp < ucomps; iComp++)
				{
					oneEntry[iComp] = 0.0;
				}
			}
			else
			{
				for (int iComp = 0; iComp < ucomps; iComp++)
				{
					oneEntry[iComp] = MeshConstants10::DRY_STATE;
				}

			}

			rv->SetTuple(idata, oneEntry);
			idata++;
		}
		delete   valBuff;
		return rv;
	}



	int num_data_layers = m_num_layers;
	if (level_center == MeshConstants10::HALF_LAYER)
	{
		num_data_layers = m_num_layers - 1;
	}
	std::string        varName(a_varName);
	int * layerStarts;
	bool is_bottom = false;
	bool is_surface = false;
	if (varName.find(m_surface_state_suffix) != string::npos)
	{
		num_data_layers = 1;
		layerStarts = new int[num_data_layers];
		layerStarts[0] = m_num_layers - 1;
		if (level_center == MeshConstants10::HALF_LAYER)
		{
			layerStarts[0]--;
		}
		is_surface = true;
	}
	else if (varName.find(m_bottom_state_suffix) != string::npos)
	{
		num_data_layers = 1;
		layerStarts = new int[num_data_layers];
		is_bottom = true;
		//this is the layer next to the bottom
		//for prism center data this is the first valid bottom prim center
		layerStarts[0] = 1;
	}
	else
	{
		layerStarts = new int[num_data_layers];
		for (int iLayers = 0; iLayers < num_data_layers; iLayers++)
		{
			layerStarts[iLayers] = iLayers;
		}
	}




	//int numDataPerLayer=m_nominal_size_per_layer[data_center];


	float * valBuff;
	int totalNumLayers = m_num_layers;

	valBuff = new float[totalNumLayers*numDataPerLayer*ncomps];


	int numOfRecord = 1;
	int nodeStart = 0;
	int timeStart = a_timeState;

	SCHISMVarPtr->set_cur(timeStart);

	if (!(SCHISMVarPtr->get(valBuff)))
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve " << a_varName << " at step " << a_timeState;
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}

	long valid_var_size = 0;
	long * node_start_index = new long[numDataPerLayer];
	long * num_data_at_layer = new long[totalNumLayers];

	for (int iLayer = 0; iLayer < m_num_layers; iLayer++)
	{
		num_data_at_layer[iLayer] = 0;
	}
	int * kbp00 = m_kbp_node[a_domainID];

	if (!(data_center.compare(MeshConstants10::EDGE)))
	{
		kbp00 = m_kbp_side[a_domainID];
	}
	else if (!(data_center.compare(MeshConstants10::ELEM)))
	{
		kbp00 = m_kbp_ele[a_domainID];
	}

	int num_layers = m_num_layers;
	if (level_center == MeshConstants10::HALF_LAYER)
	{
		num_layers--;
	}
	for (int iNode = 0; iNode < numDataPerLayer; iNode++)
	{
		node_start_index[iNode] = valid_var_size;
		valid_var_size += ((num_layers - std::max(1, kbp00[iNode]) + 1)*ncomps);
		for (int iLayer = 0; iLayer < totalNumLayers; iLayer++)
		{
			if (iLayer >= (std::max(1, kbp00[iNode]) - 1))
			{
				num_data_at_layer[iLayer]++;
			}
		}
	}

	debug1 << "size of vector data " << valid_var_size << "\n";


	vtkFloatArray *rv = vtkFloatArray::New();
	rv->SetNumberOfComponents(ucomps);
	idata = 0;

	if (varName.find(m_depth_average_suffix) == string::npos)
	{
		idata = 0;
		//total number of data =layers*nodes for a time step
		int ntuples = 0;
		for (int iLayer = 0; iLayer < num_data_layers; iLayer++)
		{
			int layer = layerStarts[iLayer];
			ntuples += num_data_at_layer[layer];
		}


		if (is_bottom || is_surface) //for surface and bottom data, same as number of 2D node
		{
			ntuples = numDataPerLayer;
		}
		rv->SetNumberOfTuples(ntuples);

		debug1 << "totla number of vect tuples " << ntuples << " " << num_data_layers << " \n";

		for (int iLayer = 0; iLayer < num_data_layers; iLayer++)
		{
			int layer = layerStarts[iLayer];

			for (int iNode = 0; iNode < numDataPerLayer; iNode++)
			{

				int bottom_layer = std::max(1, kbp00[iNode]) - 1;
				if (is_bottom)
				{
					layer = bottom_layer + 1;
					if (layer > (num_layers - 1))
					{
						layer = num_layers - 1;
					}
				}

				int start_index = node_start_index[iNode];

				if (layer >= bottom_layer)
				{
					//if(!drywet[iNode])
					if ((!drywet[iNode]) || (drywet[iNode] && (!(m_dry_wet_flag))))
					{

						for (int iComp = 0; iComp < ncomps; iComp++)
						{
							oneEntry[iComp] = valBuff[start_index + (layer + 1 - std::max(1, kbp00[iNode]))*ncomps + iComp];

						}

						for (int iComp = ncomps; iComp < ucomps; iComp++)
						{
							oneEntry[iComp] = 0.0;
						}
					}
					else
					{
						for (int iComp = 0; iComp < ucomps; iComp++)
						{
							oneEntry[iComp] = MeshConstants10::DRY_STATE;
						}
					}

					rv->SetTuple(idata, oneEntry);
					idata++;
				}
			}

		}
	}
	else
	{
		idata = 0;
		//total number of data =nodes for a time step
		int ntuples = numDataPerLayer;
		rv->SetNumberOfTuples(ntuples);

		float * averageState = new float[numDataPerLayer*ncomps];

		vectorDepthAverage(averageState, valBuff, node_start_index, a_timeState, a_domainID, ncomps, data_center, level_center);

		idata = 0;
		for (int iNode = 0; iNode < numDataPerLayer; iNode++)
		{

			//if(!drywet[iNode])
			if ((!drywet[iNode]) || (drywet[iNode] && (!(m_dry_wet_flag))))
			{
				for (int iComp = 0; iComp < ncomps; iComp++)
				{
					oneEntry[iComp] = averageState[iNode*ncomps + iComp];
				}
				for (int iComp = ncomps; iComp < ucomps; iComp++)
				{
					oneEntry[iComp] = 0.0;
				}
			}
			else
			{
				for (int iComp = 0; iComp < ucomps; iComp++)
				{
					oneEntry[iComp] = MeshConstants10::DRY_STATE;
				}
			}

			rv->SetTuple(idata, oneEntry);
			idata++;
		}
		delete averageState;
	}


	delete valBuff;
	delete layerStarts;
	delete oneEntry;
	return rv;
}

void   avtMDSCHISMFileFormatImpl::prepare_average(int * & a_kbp,
	int * & a_mapper,
	float * & a_zPtr,
	const int & a_timeState,
	const int & a_domainID,
	const std::string & a_horizontal_center,
	const std::string & a_vertical_center)
{

	a_kbp = m_kbp_node[a_domainID];

	if (a_horizontal_center == MeshConstants10::EDGE)
	{
		a_kbp = m_kbp_side[a_domainID];
		a_mapper = new int[m_num_mesh_edges[a_domainID]*m_num_layers];
		for (int i = 0; i < m_num_mesh_edges[a_domainID]*m_num_layers; i++)
		{
			a_mapper[i] = MeshConstants10::INVALID_NUM;
		}
		a_zPtr = new float[m_total_valid_3D_side[a_domainID]];
		//loadAndCacheZSide(a_timeState,a_zPtr);
		if (!m_external_mesh_providers[a_domainID]->zSideCenter3D(a_zPtr, a_timeState))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve side z at step " << a_timeState;
			EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
		}

		debug1 << " mapping edge 150 to 3d id:\n";
		int Index = 0;
		for (int iLayer = 0; iLayer < m_num_layers; iLayer++)
		{
			for (int iEdge = 0; iEdge < m_num_mesh_edges[a_domainID]; iEdge++)
			{
				int bottomLayer = a_kbp[iEdge];
				if (bottomLayer <= (iLayer + 1))
				{
					a_mapper[iLayer + iEdge * m_num_layers] = Index;
					if (iEdge == 150)
					{
						debug1 << iLayer << "->" << Index << " " << a_zPtr[Index] << " ";
					}
					Index++;
				}
			}
		}
		debug1 << "\n";
	}
	else if (a_horizontal_center == MeshConstants10::ELEM)
	{
		a_kbp = m_kbp_ele[a_domainID];
		a_mapper = new int[m_num_mesh_faces[a_domainID]*m_num_layers];
		for (int i = 0; i < m_num_mesh_faces[a_domainID]*m_num_layers; i++)
		{
			a_mapper[i] = MeshConstants10::INVALID_NUM;
		}

		int Index = 0;
		for (int iLayer = 0; iLayer < m_num_layers; iLayer++)
		{
			for (int iEle = 0; iEle < m_num_mesh_faces[a_domainID]; iEle++)
			{
				int bottomLayer = std::max(1, a_kbp[iEle]);
				if (bottomLayer <= (iLayer + 1))
				{
					a_mapper[iLayer + iEle * m_num_layers] = Index;
					Index++;
				}
			}
		}
		a_zPtr = new float[m_total_valid_3D_ele[a_domainID]];
		if (!m_external_mesh_providers[a_domainID]->zEleCenter3D(a_zPtr, a_timeState))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve ele z at step " << a_timeState;
			EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
		}
	}
	else
	{
		a_mapper = new int[m_num_mesh_nodes[a_domainID]*m_num_layers];
		for (int i = 0; i < m_num_mesh_nodes[a_domainID]*m_num_layers; i++)
		{
			a_mapper[i] = MeshConstants10::INVALID_NUM;
		}
		int Index = 0;
		for (int iLayer = 0; iLayer < m_num_layers; iLayer++)
		{
			for (int iNode = 0; iNode < m_num_mesh_nodes[a_domainID]; iNode++)
			{
				int bottomLayer = a_kbp[iNode];
				if (bottomLayer <= (iLayer + 1))
				{
					a_mapper[iLayer + iNode * m_num_layers] = Index;
					Index++;
				}
			}
		}
		a_zPtr = new float[m_total_valid_3D_point[a_domainID]];
		debug1 << "load z for average size " << m_total_valid_3D_point[a_domainID] << "\n";
		//loadAndCacheZ(a_timeState,a_zPtr);
		if (!m_external_mesh_providers[a_domainID]->zcoords3D(a_zPtr, a_timeState))
		{
			stringstream msgStream(stringstream::out);
			msgStream << "Fail to retrieve node z at step " << a_timeState;
			EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
		}
	}


}

// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::depthAverage
//
//  Purpose:
//      Averages layred state over depth and return the result
//      
//      
//
//  Arguments:
//      a_averageState   result average state
//      a_layeredState   state to be averaged
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Mon Sep  14:15:00 PDT 2012
//
// ****************************************************************************
void
avtMDSCHISMFileFormatImpl::depthAverage(float             *  a_averageState,
	float             *  a_layeredState,
	long              *  a_nodeDataStart,
	const int         &  a_timeState,
	const int         &  a_domain,
	const std::string & a_horizontal_center,
	const std::string & a_vertical_center
)
{

	long numData = m_num_mesh_nodes[a_domain];
	if (a_horizontal_center == MeshConstants10::ELEM)
	{
		numData = m_num_mesh_faces[a_domain];
	}
	else if (a_horizontal_center == MeshConstants10::EDGE)
	{
		numData = m_num_mesh_edges[a_domain];
	}
	else
	{
	}

	int * kbp;
	int * mapper;
	float * zPtr;

	prepare_average(kbp, mapper, zPtr, a_timeState, a_domain,a_horizontal_center, a_vertical_center);

	debug1 << "aveages states output \n";

	for (int iNode = 0; iNode < numData; iNode++)
	{
		if (kbp[iNode] <= m_num_layers)
		{
			// this is number of valid data
			int num_valid_data_layers = m_num_layers - std::max(1, kbp[iNode]) + 1;


			// native prism center layer 1 is duplicate of layer2
			// it still have the same number of layers as the node center data
			if (a_vertical_center == MeshConstants10::HALF_LAYER)
			{
				num_valid_data_layers--;
			}

			if (num_valid_data_layers < 1)
			{
				stringstream msgStream(stringstream::out);
				msgStream << "less than one layer at node" << iNode << " when averaging\n";
				EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
			}


			float*  states = new float[num_valid_data_layers];
			int     startid = a_nodeDataStart[iNode];

			for (int iLayer = 0; iLayer < num_valid_data_layers; iLayer++)
			{
				states[iLayer] = a_layeredState[startid + iLayer];
			}


			float*  zCoords = new float[m_num_layers - std::max(1, kbp[iNode]) + 1];  //take all the valid zs for this node

			int iiLayer = 0;
			for (int iLayer = 0; iLayer < m_num_layers; iLayer++)
			{
				if (iLayer >= (std::max(1, kbp[iNode]) - 1))
				{
					int index = mapper[iLayer + iNode * m_num_layers];
					zCoords[iiLayer] = zPtr[index];
					iiLayer++;
				}

			}



			int num_valid_z_layers = iiLayer;

			//fix for possible inactive z layers filled with 0.0 elevation 
			if (kbp[iNode] == 0)
			{
				int last_level_be_removed = 0;
				for (int ilevel = 0; ilevel < num_valid_z_layers - 1; ilevel++)
				{
					if (zCoords[ilevel + 1] < zCoords[ilevel])
					{
						last_level_be_removed = ilevel;
						//debug1<<iNode<<" has done inactive z removal to "<<last_level_be_removed<<"\n";
						break;
					}
				}
				num_valid_z_layers = num_valid_z_layers - last_level_be_removed - 1;
				num_valid_data_layers = num_valid_z_layers;
				for (int ilevel = 0; ilevel < num_valid_z_layers; ilevel++)
				{
					states[ilevel] = states[ilevel + last_level_be_removed + 1];
					zCoords[ilevel] = zCoords[ilevel + last_level_be_removed + 1];
				}

			}

			float averageState;
			if (a_vertical_center != MeshConstants10::HALF_LAYER)
			{

				averageState = trapezoidAverage(states,
					zCoords,
					num_valid_data_layers);
			}
			else
			{
				averageState = rectAverage(states,
					zCoords,
					num_valid_z_layers);
			}

			a_averageState[iNode] = averageState;

			if (iNode == 20002)
			{
				debug1 << "var vals for " << iNode << " averageing:\n";
				debug1 << "level center " << a_vertical_center << "\n";
				debug1 << " bottom level:" << kbp[iNode] << "\n";
				debug1 << " valid num of level:" << num_valid_data_layers << "\n";
				for (int iLayer = 0; iLayer < num_valid_data_layers; iLayer++)
				{
					debug1 << states[iLayer] << " " << zCoords[iLayer] << "\n";
				}
				debug1 << "averaged to " << averageState << "\n";
			}
			delete   states;
			delete   zCoords;
		}
		else
		{
			a_averageState[iNode] = MeshConstants10::DRY_STATE;
		}
	}
	delete zPtr;
	delete mapper;
	//debug1<<numData<<" center "<<m_data_center<<"\n";
}


// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::vectorDepthAverage
//
//  Purpose:
//      Averages layred state over depth and return the result
//      
//      
//
//  Arguments:
//      a_averageState   result average state
//      a_layeredState   state to be averaged
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Mon Sep  14:15:00 PDT 2012
//
// ****************************************************************************
void
avtMDSCHISMFileFormatImpl::vectorDepthAverage(float        *  a_averageState,
	float        *  a_layeredState,
	long          *  a_nodeDataStart,
	const int    &  a_timeState,
	const int    &  a_domain,
	const int    &  a_ncomps,
	const std::string & a_horizontal_center,
	const std::string & a_vertical_center
)
{

	long numData = m_num_mesh_nodes[a_domain];
	if (a_horizontal_center == MeshConstants10::ELEM)
	{
		numData = m_num_mesh_faces[a_domain];
	}
	else if (a_horizontal_center == MeshConstants10::EDGE)
	{
		numData = m_num_mesh_edges[a_domain];
	}
	else
	{
	}

	int * kbp;
	int * mapper;
	float * zPtr;

	prepare_average(kbp, mapper, zPtr, a_timeState, a_domain,a_horizontal_center, a_vertical_center);

	for (int iNode = 0; iNode < numData; iNode++)
	{
		if (kbp[iNode] <= m_num_layers)
		{
			int num_valid_data_layers = m_num_layers - std::max(1, kbp[iNode]) + 1;

			if (a_vertical_center == MeshConstants10::HALF_LAYER)
			{
				num_valid_data_layers--;
			}

			if (num_valid_data_layers < 1)
			{
				stringstream msgStream(stringstream::out);
				msgStream << "less than one level at node" << iNode << " when averaging\n";
				EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
			}


			float**  states = new float*[a_ncomps];
			for (int iComp = 0; iComp < a_ncomps; iComp++)
			{
				states[iComp] = new float[num_valid_data_layers];
			}
			int     startid = a_nodeDataStart[iNode];
			for (int iLayer = 0; iLayer < num_valid_data_layers; iLayer++)
			{
				for (int iComp = 0; iComp < a_ncomps; iComp++)
				{
					states[iComp][iLayer] = a_layeredState[iComp + iLayer * a_ncomps + startid];
					if (iNode == 1850)
					{
						debug1 << states[iComp][iLayer] << " ";
					}
				}
			}

			float*  zCoords = new float[m_num_layers - std::max(1, kbp[iNode]) + 1];  //take all the zs for the mesh  

			int iiLayer = 0;

			if (iNode == 1850)
			{
				debug1 << " z of 1850 : ";
			}


			for (int iLayer = 0; iLayer < m_num_layers; iLayer++)
			{
				if (iLayer >= (std::max(1, kbp[iNode]) - 1))
				{
					long index = mapper[iLayer + iNode * m_num_layers];
					zCoords[iiLayer] = zPtr[index];
					if (iNode == 1850)
					{
						debug1 << zCoords[iiLayer] << " " << iLayer << "->" << index << " ";
					}
					iiLayer++;
				}

			}

			int num_valid_z_layers = iiLayer;
			//fix for possible inactive z layers filled with 0.0 elevation 
			if (kbp[iNode] == 0)
			{
				int last_level_be_removed = 0;
				for (int ilevel = 0; ilevel < num_valid_z_layers - 1; ilevel++)
				{
					if (zCoords[ilevel + 1] < zCoords[ilevel])
					{
						last_level_be_removed = ilevel;
						break;
					}
				}
				num_valid_z_layers = num_valid_z_layers - last_level_be_removed - 1;
				num_valid_data_layers = num_valid_z_layers;
				for (int ilevel = 0; ilevel < num_valid_z_layers; ilevel++)
				{
					for (int iComp = 0; iComp < a_ncomps; iComp++)
					{
						states[iComp][ilevel] = states[iComp][ilevel + last_level_be_removed + 1];
					}
					zCoords[ilevel] = zCoords[ilevel + last_level_be_removed + 1];
				}

			}

			if (iNode == 1850)
			{
				debug1 << " averaged sate of 1850 : ";
			}


			for (int iComp = 0; iComp < a_ncomps; iComp++)
			{
				float averageState;
				if (a_vertical_center != MeshConstants10::HALF_LAYER)
				{
					averageState = trapezoidAverage(states[iComp],
						zCoords,
						num_valid_data_layers);


				}
				else
				{
					averageState = rectAverage(states[iComp],
						zCoords,
						num_valid_z_layers);
					if (iNode == 1850)
					{
						debug1 << averageState << " ";
					}

				}

				a_averageState[iNode*a_ncomps + iComp] = averageState;
			}
			if (iNode == 1850)
			{
				debug1 << " end 1850\n ";
			}
			for (int iComp = 0; iComp < a_ncomps; iComp++)
			{
				delete  states[iComp];
			}
			delete   states;
			delete   zCoords;
		}
		else
		{
			for (int iComp = 0; iComp < a_ncomps; iComp++)
			{
				a_averageState[iNode*a_ncomps + iComp] = MeshConstants10::DRY_STATE;
			}
		}
	}

	delete zPtr;
	delete mapper;

}


// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::GetSingleLayerVar
//
//  Purpose:
//      Gets a scalar variable asscocated with mesh with only one layer.
//      
//      
//
//  Arguments:
//      a_timeState The index of the a_timeState.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      a_varName   The netcdf name of the variable requested.
//
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Wed Aug 21 10:11:04 PDT 2012
//
// ****************************************************************************

void
avtMDSCHISMFileFormatImpl::getSingleLayerVar(float    *          a_valBuff,
	SCHISMFile10*         a_SCHISMOutPtr,
	const int &         a_timeState,
	const std::string&  a_varName) const
{
	time_t startTicks = clock();
	SCHISMVar10 * SCHISMVarPtr = a_SCHISMOutPtr->get_var(a_varName);
	if (!(SCHISMVarPtr->is_valid()))
	{

		EXCEPTION1(InvalidVariableException, a_varName);
	}
	debug1 << "begin to read " << a_varName << endl;

	int nodeIndex0 = 0;
	int timeRecord = a_timeState;

	SCHISMVarPtr->set_cur(timeRecord);

	debug1 << "set start of  " << a_varName << endl;
	int numOfRecord = 1;

	if (!SCHISMVarPtr->get(a_valBuff))
	{
		stringstream msgStream(stringstream::out);
		msgStream << "Fail to retrieve " << a_varName << " at step " << a_timeState;
		EXCEPTION3(DBYieldedNoDataException, m_data_file, m_plugin_name, msgStream.str());
	}

	time_t endTicks = clock();
	debug1 << "time used in getting var " << a_varName << ":" << endTicks - startTicks << endl;
}




// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::ActivateTimestep
//
//  Purpose:
//      
//      
//  Programmer: qshu -- generated by xml2avt
//  Creation:   Mon Aug 6 09:53:36 PDT 2012
//
// ****************************************************************************

//void  avtMDSCHISMFileFormatImpl::ActivateTimestep(const std::string& a_filename)
//{
//	Initialize(a_filename);
//}

void    avtMDSCHISMFileFormatImpl::count_node_side_num_domain(const std::string& a_path,const int& num_proc)
{
	//m_side_num_domain = new int[m_global_num_side];
	//m_node_num_domain = new int[m_global_num_node];

	std::string local_global_file;
    //std:string pt = "\\";


	int myrank = 0;
#ifdef PARALLEL
	myrank = PAR_Rank();
#endif

	//if (myrank == 0)
	{

		for (int i = 0; i < num_proc; i++)
		{
			std::stringstream ss;
			ss << std::setw(4) << std::setfill('0') << i;
			std::string domain_str = ss.str();
			
#ifdef _WIN32
			local_global_file = a_path + "\\local_to_global_"+ domain_str;
#else
			local_global_file = a_path + "/local_to_global_"+ domain_str;
#endif
			ifstream*    localFileStream = new ifstream(local_global_file.c_str());
			debug1 << "reading global id from " << local_global_file << "\n";
			if (!localFileStream->good())
			{
				EXCEPTION1(InvalidFilesException,local_global_file+" is not valid");
			}
			std::string  lineTemp;

			std::getline(*localFileStream, lineTemp);
			std::getline(*localFileStream, lineTemp);
			long num_element;
			*(localFileStream) >> num_element;
			
			for (long iEle = 0; iEle < num_element; iEle++)
			{
				long v1, v2;
				*(localFileStream) >> v1 >> v2;
				
			}
			long num_node;
			*(localFileStream) >> num_node;
			

			for (long iNode = 0; iNode < num_node; iNode++)
			{
				long v1, v2;
				*(localFileStream) >> v1 >> v2;
				m_node_num_domain[v2 - 1]++;
				
			}

			long num_side;
			*(localFileStream) >> num_side;
			
			for (long iSide = 0; iSide < num_side; iSide++)
			{
				long v1, v2;
				*(localFileStream) >> v1 >> v2;
				m_side_num_domain[v2 - 1]++;
			}
			localFileStream->close();
			delete localFileStream;
		}

	}

}

int   avtMDSCHISMFileFormatImpl::load_per_proc_file(const std::string& a_path, int & num_node, int & num_side, int & num_ele) const
{

	std::string local_global_file;
#ifdef _WIN32
	local_global_file = a_path + "\\local_to_global_0000";
#else
	local_global_file = a_path + "/local_to_global_0000";
#endif


	ifstream*    localFileStream = new ifstream(local_global_file.c_str());
	if (!localFileStream->good())
	{
	
	   EXCEPTION1(InvalidDBTypeException, "not valid schsim NC output");
	
	}
	std::string  lineTemp;

	std::getline(*localFileStream, lineTemp);
    std:stringstream lstream(lineTemp);
	
	int nl, nproc;
	lstream >> num_side >> num_ele >> num_node >> nl >> nproc;
	localFileStream->close();
	delete localFileStream;
	return nproc;

}
int avtMDSCHISMFileFormatImpl::num_domain()
{

	return m_number_domain;
}
// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::Initialize
//
//  Purpose:
//      Try to open the saved input data filename using ncfile and check
//      if the file is a valid netcdf file. If valid save ncfile pointer.
//     
//
//  Programmer: qshu
//  Creation:   Mon Aug 6 09:53:36 PDT 2012
//
// ****************************************************************************
void avtMDSCHISMFileFormatImpl::Initialize(std::string a_data_file)
{

	if (!m_initialized)
	{
		bool okay = false;
		// Open the file specified by the filename argument here using
		// ncfile API. See if the file has the right things in
		// it. If so, set okay to true.
		debug1 << "file is going to be opened "<<m_data_file<<"\n";
		m_data_file = a_data_file;
		
		size_t found = m_data_file.find_last_of("/\\");
		debug1 << "found is  " << found<< "\n";
		size_t found1 = m_data_file.find_first_of("_");
		size_t found2 = m_data_file.find_last_of("_");
		
		m_data_file_path = m_data_file.substr(0, found);
		
		//std::string data_file_part1 = m_data_file.substr(0, found1);
		std::string data_file_part2 = m_data_file.substr(found2);
		//debug1 << "path:" << m_data_file_path << "\n";
		debug1 << "part2:" << data_file_part2 << "\n";
		//debug1 << m_data_file_path + "\\" + "schout_" + "0000" + "_" + data_file_part2+"\n";
		int nproc = 1;
		int myrank = 0;
#ifdef PARALLEL
		nproc = PAR_Size();
	    myrank = PAR_Rank();
#endif
		int ndomain = 0;
		if (myrank == 0)
		{
			ndomain = load_per_proc_file(m_data_file_path, m_global_num_node, m_global_num_side, m_global_num_ele);
			
		}
		debug1 << "global num before broadcast:" << m_global_num_node << " " << m_global_num_ele << " " << m_global_num_side << "\n";
#ifdef PARALLEL
		BroadcastInt(ndomain);
		BroadcastInt(m_global_num_node);
		BroadcastInt(m_global_num_ele);
		BroadcastInt(m_global_num_side);
		debug1 << "global num after broadcast:" << m_global_num_node << " " << m_global_num_ele << " " << m_global_num_side << "\n";
#endif

		for (int i = 0; i < m_global_num_side; i++)
		{
			m_side_num_domain.push_back(0);
		}
		for (int i = 0; i < m_global_num_node; i++)
		{
			m_node_num_domain.push_back(0);
		}

	
		if (myrank == 0)
		{
			count_node_side_num_domain(m_data_file_path, ndomain);
		}
#ifdef PARALLEL

		BroadcastIntVector(m_side_num_domain, 0);
		BroadcastIntVector(m_node_num_domain, 0);

#endif

#ifdef PARALLEL

		//int *buff1 = new int[20000];
		//int *buff2 = new int[20000];
		//if (myrank == 0)
		//{
		//	for (int i = 0; i < 20000; i++)
		//		buff1[i] = m_side_num_domain[i];
		//	for (int i = 0; i < 20000; i++)
		//		buff2[i] = m_node_num_domain[i];
		//	debug1 << "broadcast domain num arr\n";
		//	int bsize1, bsize2;
		//	bsize1 = m_global_num_node;
		//	bsize2 = m_global_num_side;
		//	debug1 << "global num after broadcast:" << bsize1<<" "<<bsize2<< "\n";
		//	BroadcastIntArray(buff1, 20000);
		//	BroadcastIntArray(buff2, 20000);
		//	debug1 << "finished broadcast domain num arr\n";
		//}
		//if (myrank != 0)
		//{
		//	for (int i = 0; i < m_global_num_side; i++)
		//		m_side_num_domain[i]=buff1[i];
		//	for (int i = 0; i < m_global_num_node; i++)
		//		m_node_num_domain[i]=buff2[i];
		//}
		//delete[] buff1;
		//delete[] buff2;
#endif PARALLEL

		/*std::string d1 = m_data_file_path + "\\schout_0000_5.nc";
		std::string l1 = m_data_file_path +  "\\local_to_global_0000";
		MDSchismOutput* a_data_file_ptr = NULL;
	    debug1 << "creating file ptr " << 0 << "\n";
	    a_data_file_ptr = new MDSchismOutput(d1, l1);*/
		m_number_domain = ndomain;
		int ndomain_per_thread = ndomain / nproc;

		int oneExtraUntil= ndomain % nproc;
	
	    std:string pt = "\\";
		debug1 << "proc file be opened ndomain:" << ndomain << "\n";
#ifdef _WIN32
		pt = "\\";
#else
		pt = "/";
#endif

#ifdef PARALLEL
		

		int num_domain_resides = ndomain_per_thread + (myrank < oneExtraUntil ? 1 : 0);
		int *all_domain_resides = new int[num_domain_resides];
		
		int start_domain = 0;
		
	
		for (int i = 0; i < ndomain_per_thread; i++)
		{
			all_domain_resides[i] = start_domain + i;
		}
		
		for (int i = 0; i < nproc; i++)
		{
			if (i == myrank)
			{
				for (int j = 0; j < num_domain_resides; j++)
				{
					all_domain_resides[j]=j + start_domain;
				}
			}
			start_domain += ndomain_per_thread + (i < oneExtraUntil ? 1 : 0);
		}

		debug1 << "mesh has been distributed num domain resides:"<<num_domain_resides<<"\n";

		for (int i = 0; i < num_domain_resides; i++)
		{
			int id = all_domain_resides[i];
			std::stringstream ss;
			ss << std::setw(4) << std::setfill('0') << id;
			std::string domain_str = ss.str();
			//std::string domain_data_file = data_file_part1 + domain_str + "_" + data_file_part2;
			std::string domain_data_file = m_data_file_path+ pt +"schout_"+ domain_str + data_file_part2;
			debug1 << "loading " << domain_data_file << "\n";
			std::string local_mesh_file = m_data_file_path + pt + "local_to_global_" + domain_str;
			MDSchismOutput* a_data_file_ptr = NULL;
			try
			{
				debug1 << "creating file ptr "<<i << "\n";
				a_data_file_ptr = new MDSchismOutput(domain_data_file, local_mesh_file);
				debug1 << "created file ptr " << a_data_file_ptr << "\n";
			}
			catch (SCHISMFileException10 e)
			{
				EXCEPTION1(InvalidDBTypeException, e.what());
			}
			catch (...)
			{
				EXCEPTION1(InvalidDBTypeException, "not valid schsim NC output");
			}
			m_data_files[id] = a_data_file_ptr;

			MDSCHISMMeshProvider * a_local_mesh_ptr = NULL;
			try
			{
				a_local_mesh_ptr = new MDSCHISMMeshProvider(domain_data_file, a_data_file_ptr, local_mesh_file);
				//a_local_mesh_ptr->set_data_file(a_data_file_ptr);
				debug1 << "created mesh ptr " << a_local_mesh_ptr << "\n";
			}
			catch (SCHISMFileException10 e)
			{
				EXCEPTION1(InvalidDBTypeException, e.what());
			}
			catch (...)
			{
				EXCEPTION1(InvalidDBTypeException, "no valid mesh exist");
			}
			m_external_mesh_providers[id] = a_local_mesh_ptr;
		}
#else
		for (int id = 0; id < ndomain; id++)
		{

			std::stringstream ss;
			ss << std::setw(4) << std::setfill('0') << id;
			std::string domain_str = ss.str();
			std::string domain_data_file = m_data_file_path + pt + "schout_" + domain_str  + data_file_part2;
			std::string local_mesh_file = m_data_file_path + pt + "local_to_global_" + domain_str;
			MDSchismOutput* a_data_file_ptr = NULL;
			try
			{
				a_data_file_ptr = new MDSchismOutput(domain_data_file, local_mesh_file);
			}
			catch (SCHISMFileException10 e)
			{
				EXCEPTION1(InvalidDBTypeException, e.what());
			}
			catch (...)
			{
				EXCEPTION1(InvalidDBTypeException, "not valid schsim NC output");
			}
			m_data_files[id] = a_data_file_ptr;

			MDSCHISMMeshProvider * a_local_mesh_ptr = NULL;
			try
			{
				a_local_mesh_ptr = new MDSCHISMMeshProvider(domain_data_file, a_data_file_ptr, local_mesh_file);
				//a_local_mesh_ptr->set_data_file(a_data_file_ptr);
			}
			catch (SCHISMFileException10 e)
			{
				EXCEPTION1(InvalidDBTypeException, e.what());
			}
			catch (...)
			{
				EXCEPTION1(InvalidDBTypeException, "no valid mesh exist");
			}
			m_external_mesh_providers[id] = a_local_mesh_ptr;
	     }
#endif


		debug1 << "file is opened\n";

		// need to add get_dry_wetl flag to MDSchismOuput class
		//m_dry_wet_flag = m_data_files.begin()->second->get_dry_wet_val_flag();
		debug1 << "wet_dry_flag is " << m_dry_wet_flag << "\n";
		debug1 << "begin get dim\n";

		getTime();
		

		debug1 << "getting dim";

		std::map<int, MDSCHISMMeshProvider *>::iterator it;

		for (it = m_external_mesh_providers.begin(); it != m_external_mesh_providers.end(); it++)
		{
			getMeshDimensions(it->first,it->second);
		}

		
		debug1 << "got dimension\n";

		debug1 << "loading coords";

		m_mesh_is_static = m_external_mesh_providers.begin()->second->mesh3d_is_static();

		PopulateVarMap();

		int current_time = 0;
		for (it = m_external_mesh_providers.begin(); it != m_external_mesh_providers.end(); it++)
		{
			this->load_ele_dry_wet(current_time,it->first);
			//debug1 << "done load dry/wet for domain:" << it->first << "\n";
			this->load_bottom(current_time,it->first);
			debug1 << "done load bottom for domain:" << it->first << "\n";
			m_kbp_prism[it->first] = new int [m_num_mesh_faces[it->first]];
		}

		m_initialized = true;

		debug1 << "done initialize\n";
	}
}




// a reduant populate var map dic, fix a bug in which this dic is 
// lost when switch a different data file group in the same .visit file
// called in initialize()
void avtMDSCHISMFileFormatImpl::PopulateVarMap()
{
	int numVar = m_data_files.begin()->second->num_vars();
	//std::string  location = m_data_file_ptr->data_center();

	for (int iVar = 0; iVar < numVar; iVar++)
	{
		SCHISMVar10*  varPtr = m_data_files.begin()->second->get_var(iVar);
		std::string varName = varPtr->name();
		
		std::string  location = varPtr->get_horizontal_center();
		//debug1 << "loading var name " << varName << "\n";
		if ((varName == m_node_surface) || (varName == m_node_depth))
		{
			continue;
		}
		std::string  label;
		label = varName;
		// this dic make it easy to find out data set for a visit plot variable
		m_var_name_label_map[label] = varName;

		if ((location == NODE) || (location == FACE))
		{
			// all those surface, botootm and average are based on original data set
			m_var_name_label_map[label + m_surface_state_suffix] = varName;
			m_var_name_label_map[label + m_bottom_state_suffix] = varName;
			m_var_name_label_map[label + m_depth_average_suffix] = varName;

			m_var_mesh_map[label + m_surface_state_suffix] = m_mesh_2d;
			m_var_mesh_map[label + m_bottom_state_suffix] = m_mesh_2d;
			m_var_mesh_map[label + m_depth_average_suffix] = m_mesh_2d;

		}
		// omit unkown center data
		else
		{
			continue;
		}
	}
	//debug1 << "done loading all var name \n";
}






// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::GetTimes
//
//  Purpose:
//      overloaded public interface to return time steps
//      
//     
//
//  Programmer: qshu
//  Creation:   Mon Aug 15 03:02:00 PDT 2012
//
// ****************************************************************************
void   avtMDSCHISMFileFormatImpl::GetTimes(std::vector<double> & a_times)
{

	//copy saved time into a_times
	for (int i = 0; i < m_num_time_step; i++)
	{
		a_times.push_back(m_time_ptr[i]);
	}

}

// ****************************************************************************
//  Method: avtMDSCHISMFileFormatImpl::getTime
//
//  Purpose:
//      Try to get time steps saved in this file
//      
//     
//
//  Programmer: qshu
//  Creation:   Mon Aug 15 10:11:00 PDT 2012
//
// ****************************************************************************
void  avtMDSCHISMFileFormatImpl::getTime()
{

	SCHISMDim10 * dimTimePtr = m_data_files.begin()->second->get_dim(m_dim_time);

	if (dimTimePtr->is_valid())
	{
		m_num_time_step = dimTimePtr->size();
	}
	else
	{
		EXCEPTION1(InvalidVariableException, m_dim_time);
	}


	if (m_num_time_step)
	{
		m_time_ptr = new float[m_num_time_step];

		SCHISMVar10* SCHISMTimePtr = m_data_files.begin()->second->get_var(m_time);

		if (SCHISMTimePtr->is_valid())
		{
			SCHISMTimePtr->get(m_time_ptr);
		}
		else
		{
			EXCEPTION1(InvalidVariableException, m_time);
		}

	}
}



bool  avtMDSCHISMFileFormatImpl::SCHISMVarIs3D(SCHISMVar10*  a_varPtr) const
{
	int totalNumDim = a_varPtr->num_dims();
	for (int dim = 0; dim < totalNumDim; dim++)
	{
		SCHISMDim10* dimPtr = a_varPtr->get_dim(dim);
		if (dimPtr->name() == m_dim_layers)
			return true;
	}
	return false;
}
bool  avtMDSCHISMFileFormatImpl::SCHISMVarIsVector(SCHISMVar10* a_varPtr) const
{
	int totalNumDim = a_varPtr->num_dims();
	for (int dim = 0; dim < totalNumDim; dim++)
	{
		SCHISMDim10* dimPtr = a_varPtr->get_dim(dim);
		if (dimPtr->name() == m_dim_var_component)
			return true;
	}
	return false;
}

static Registrar registrar("md_nc4", &avtMDSCHISMFileFormatImpl::create);
